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Abstract 

Data  Assimilation  for  Ocean  Acoustic  Tomography 

by  Chris  G.  Walter 

Chair  of  Supervisory  Committee 

Research  Associate  Professor  James  A.  Mercer 
Department  of  Geophysics 

This  research  evaluates  the  effectiveness  of  assimilating  ocean  acoustic  tomog¬ 
raphy  data  into  numerical  ocean  models.  Acoustic  tomography  uses  sound  to 
remotely  sample  integrated  or  averaged  properties  of  the  ocean.  These  prop¬ 
erties  include  soundspeed  and  current  velocity.  Three  different  models  are 
used  in  the  assimilation:  1)  a  linear  Rossby  wave  model  without  advection, 
2)  a  Rossby  wave  model  with  advection,  and  3)  a  non-linear  quasi-geostrophic 
model.  We  examine  the  tomographic  data’s  effectiveness  by  assimilating  non¬ 
averaging  point  measurements,  such  as  temperature  data,  and  integral  data 
with  the  Rossby  wave  model.  First,  simulated  data  are  used  in  benchmark 
experiments.  These  are  followed  by  a  second  series  of  experiments  in  which 
actual  tomographic  data  collected  during  the  Acoustic  Mid-Ocean  Dynamics 
Experiment  (AMODE)  are  used.  We  compare  the  value  of  the  integral  mea¬ 
surements  in  terms  of  forecasting  accuracy  in  both  physical  and  spectral  space. 
Although  both  forms  of  data  constrain  the  model  equally  as  well  in  an  average 
sense  (the  number  of  data  were  chosen  to  do  this),  the  tomographic  data  are 
found  to  constrain  low  wavenumber  components  of  the  model  more  accurately. 
Differences  between  the  model  estimates  and  data  are  within  expected  error 
bars.  These  error  bars  are  based  on  prior  statistical  assumptions  about  errors 
in  the  model  and  observations,  boundary  conditions  and  initial  conditions.  The 
Rossby  wave  models  with  and  without  advection  account  for  62.1%  and  60.7% 
of  the  data  variance  respectively  and  the  quasi-geostrophic  model  accounts  for 


66.0%  of  the  variance.  The  purpose  of  this  research  is  not  to  evaluate  any 
particular  numerical  model,  but  rather  to  better  understand  the  usefulness  of 
assimilating  tomographic  measurements  into  numerical  models. 


FOREWORD 


This  report  is  a  slightly  revised  version  of  the  dissertation  submitted  in  partial  fulfillment 
of  the  requirements  for  the  Doctor  of  Philosophy  degree  at  the  University  of  Washington  in  June 
1999.  Dr.  James  A.  Mercer  was  the  chair  of  the  supervisory  committee. 
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Chapter  1 

INTRODUCTION 


1.1  Motivation 

Measuring  and  predicting  the  ocean’s  mesoscale  circulation  presents  a 
formidable  task  for  oceanographers.  Ocean  mesoscale  features  are  compara¬ 
ble  to  storms  and  weather  fronts  in  the  atmosphere  and  are  of  great  interest 
to  a  wide  variety  of  disciplines,  ranging  from  global  meteorology  to  underwater 
tracking  and  submarine  detection.  In  order  to  measure  and  adequately  resolve 
these  mesoscale  features,  their  length  scales  must  be  monitored  on  the  order 
of  the  Rossby  radius  («  50  km  in  the  mid-latitude  North  Atlantic)  and  time 
periods  on  the  order  of  weeks  (Hogg,  1996).  This  places  severe  demands  on 
most  ocean  measurement  systems.  One  system  that  is  particularly  efficient  at 
measuring  features  of  this  space  and  time  scale  is  ocean  acoustic  tomography 
(Munk  and  Wunsch,  1979). 

With  ocean  tomography,  acoustic  pulses  are  transmitted  through  the  ocean’s 
interior  using  an  array  of  sources  and  receivers.  The  tomographic  data  consist 
of  travel  times  or  “times  of  flight”.  From  these  travel  times,  it  is  possible  to 
deduce  the  horizontally  and  vertically  averaged  soundspeed  and  currents  of 
the  intervening  ocean.  This  averaging  makes  the  measurements  integral  in 
nature,  and  because  of  this  integral  nature,  the  data  are  more  sensitive  to  eddy 
and  mesoscale  features. 

Typical  non-tomographic  measurements  consist  of  in-situ  observations 
such  as  CTD  (Conductivity  Temperature  and  Depth)  and  XBT  (expendable 
BathyThermograph)  casts.  These  forms  of  data  as  well  as  satellite  altime¬ 
try  data  can  be  classified  as  “point”  or  non-averaging  measurements.  (The 
satellite  data  actually  measure  vertical  but  not  horizontal  averages  of  the  wa- 
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ter  density).  Point  measurements,  such  as  these,  are  sensitive  to  all  horizon¬ 
tal  length  scales  (Cornuelle,  1996).  This  broad  sensitivity  makes  them  prone 
to  contamination  from  short  length  scale  noise  processes  (such  as  internal  or 
gravity  waves).  Figure  1.1  helps  to  clarify  these  concepts.  It  shows  various 
measurement  techniques  and  the  domains  over  which  the  measurements  inte¬ 
grate  or  average  during  the  sampling  process  (aside  from  any  post-processing 
or  smoothing).  It  is  important  to  note  that  noise  from  these  non-integrated 
domains  can  contaminate  the  measurements. 

In  order  to  best  estimate  the  ocean  properties,  data  assimilation  or  inverse 
methods  are  used  to  reconstruct  the  temperature  (soundspeed)  and  current 
fields  from  the  data.  The  goal  in  the  assimilation  is  to  take  travel  time  ob¬ 
servations  that  are  sparse  and  possibly  contaminated  with  noise,  and  produce 
the  best  estimate  of  the  ocean’s  4-dimensional  circulation.  Data  assimilation  is 
particularly  important  with  data  types  such  as  tomography.  The  non-local  na¬ 
ture  of  tomography  can  result  in  significant  cross  correlations  in  the  properties 
that  tomography  is  trying  to  measure  (Cornuelle,  1996).  If  these  cross  correla¬ 
tions  are  not  properly  taken  into  account,  information  will  be  lost  or  neglected 
in  the  model’s  final  estimates.  In  the  past,  we  simplified  this  inverse  problem 
by  either  neglecting  the  time  dependence  of  the  data  altogether  or  by  applying 
some  simplified  time-lag  filter  to  the  data.  The  interested  reader  is  directed  to 
Munk  and  Wunsch  (1982)  or  Munk  and  Wunsch  (1979)  or  Howe  et  al.  (1987). 
Neither  of  these  approaches  makes  use  of  a  priori  knowledge  of  ocean  dynamics 
to  constrain  how  the  temperature  and  current  fields  evolve  and  may  produce 
unrealistic  estimates  of  the  circulation. 

To  help  evaluate  the  success  of  our  assimilation  methods,  we  make  use  of 
a  number  of  Twin  Experiments  (TEs).  The  two  major  portions  of  a  TE  are  a 
control  model  and  an  assimilation  estimate  of  the  control  model.  The  assimi¬ 
lation  estimate  uses  data  simulated  from  the  control  ocean.  The  control  model 
can  be  thought  of  as  the  “true  ocean”,  and  the  data  as  “true”  data,  which  may 
be  artificially  contaminated  with  noise  before  they  are  used  in  the  assimilation 
process.  Comparisons  are  made  between  the  assimilation  ocean  and  the  con¬ 
trol  ocean  to  measure  goodness  of  the  forecasts.  See  Figure  1.2  which  shows  a 
schematic  diagram  of  a  typical  twin  experiment. 
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In  Figure  1.2  we  show  the  four  major  components  that  any  ocean  data  as¬ 
similation  experiment,  whether  it  uses  simulated  data  or  real  data,  consists 
of:  1)  A  dynamical  model  which  can  evolve  an  ocean  field  in  time.  2)  Obser¬ 
vations  or  data  with  measurement  error,  often  irregularly  sampled  in  space 
and  time.  3)  Error  statistics  for  both  the  model  and  the  data.  4)  A  system  for 
combining  model  output  with  data.  The  dynamical  model  is  a  form,  almost  al¬ 
ways  greatly  simplified,  of  the  equations  of  motion  with  the  necessary  initial 
and  boundary  conditions.  In  addition  to  boundary  conditions,  measurements  in 
the  interior  (e.g.,  current  meter  data,  CTD  profiles,  acoustic  travel  times,  and 
drifter  positions)  are  generally  used.  There  are  two  fundamental  data  types: 
point  data  and  integral  data.  Examples  of  the  first  include  CTD  profiles,  pres¬ 
sure  data,  and  current  meter  data.  Examples  of  integral  measurements  are 
acoustic  travel  time  (inverted  echo  sounders  as  well  as  "tomographic"  travel 
times),  electric  field  measurements,  drifters  and  floats,  and  altimetric  sea  sur¬ 
face  height  measurements  (height  and  pressure  are  functions  of  the  integral 
of  the  water  density).  Float  or  drifter  data  are  complicated  because  they  are 
an  integral  function  of  time  and  spatially  varying  velocity,  since  they  measure 
fto  u(x(t))dt  . 

A  prescription  of  the  error  statistics  is  absolutely  necessary  for  rigorous  data 
assimilation  (Comuelle,  1996).  The  data  errors  may  be  correlated  with  one  an¬ 
other  (such  as  altimeter  orbit  errors  or  XBT  fall  rate  errors),  with  the  model 
errors,  or  with  the  model  state  itself  and  may  be  difficult  to  estimate.  Model 
errors  are  especially  difficult  to  prescribe,  because  doing  so  usually  implies  one 
has  a  better  model  available.  The  specification  of  errors  is  often  the  most  diffi¬ 
cult  (and  often  subjective)  part  of  data  assimilation. 

1.2  Background 

Acoustic  tomography  was  first  seriously  discussed  by  Munk  and  Wunsch  (1979, 
1982).  Munk  and  Wunsch  discuss  the  wisdom  of  using  information  from  acous¬ 
tic  transmissions  as  a  means  of  remotely  measuring  the  ocean  mesoscale  struc¬ 
ture.  These  acoustic  transmissions  measure  certain  properties  of  the  ocean  us¬ 
ing  travel  times  (times  of  flight)  between  distant  moorings.  By  using  an  array  of 
moorings  (actually  acoustic  transceivers),  synoptic  collections  of  “tomographic” 
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observations  may  be  made  over  vast  regions  of  the  ocean. 

Within  the  last  fifteen  years  numerous  tomographic  experiments  have  been 
conducted  in  many  oceans  of  the  world.  Excellent  reviews  of  these  experiments 
are  given  by  Munk  et  al.  (1995)  and  also  by  Worcester  et  al.  (1991).  In  ad¬ 
dition,  Dushaw  et  al.  (1993)  offer  a  valuable  summary  of  tomographic  use  of 
reciprocal  travel  time  transmissions.  Even  with  this  proliferation  of  documen¬ 
tation  extolling  tomography’s  value,  the  oceanography  community  has  shown 
reluctance  in  utilizing  the  data.  Partly,  the  reluctance  stems  from  the  dras¬ 
tic  difference  between  point  and  tomographic  data.  At  this  time,  the  intrinsic 
value  of  the  tomographic  data  is  not  understood  or  appreciated.  Data  assimi¬ 
lation  may  help  bring  tomographic  data  into  a  more  acceptable  light.  This  will 
be  achieved  by  using  data  assimilation  as  a  tool  to  learn  how  the  line-integral 
nature  of  tomographic  data  constrains  the  circulation  model.  A  major  goal  of 
this  research  is  to  demonstrate  the  complementary  nature  of  tomographic  and 
point  measurements.  They  should  be  individual  pieces  of  a  total  ocean  observ¬ 
ing  system. 

Howe  et  al.  (1987)  describe  a  method  used  to  invert  travel  times  for  sound- 
speed  and  current.  These  sum  and  difference  travel  times  measure  integral 
properties  of  soundspeed  (temperature)  and  current.  Tomography’s  integrating 
nature  acts  as  a  spatial  low  pass  filter,  and  is  more  sensitive  to  low  wavenumber 
features.  Comuelle  et  al.  (1989)  and  Comuelle  and  Howe  (1987)  have  examined 
these  measurement  properties  for  time  independent  mapping. 

The  most  common  data  analysis  methods  that  the  tomographic  community 
has  used  to  date  employ  simple  linear  inverse  methods.  These  inverse  methods 
assume  that  all  travel  time  data  is  collected  over  a  fixed  time  window.  Several 
tomographic  researchers  have  tried  limited  variations  on  this  theme.  For  ex¬ 
ample,  Howe  et  al.  (1987)  applied  a  simple  time-lag  filter  using  travel  times 
which  were  collected  with  a  single  source-receiver  pair.  Spiesberger  Metzger 
(1989)  analyzed  data  by  using  a  generalized  mapping  method  with  a  time  and 
space  varying  covariance  function.  Neither  method  considered  any  knowledge 
of  ocean  circulation  dynamics  to  constrain  the  temperature  or  current  fields. 
Several  analysis  methods,  e.g.,  those  by  Chiu  and  Desaubies  (1987)  and  Cor- 
nuelle  et  al.  (1985)  use  simple  linear  dynamics.  These  dynamics  assume  the 
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model  perturbations  are  linear  Rossby  waves  with  time-independent  ampli¬ 
tudes. 

The  term  “data  assimilation”  stems  from  the  meteorological  field.  About 
40  years  ago,  the  phrase  was  coined  as  a  method  of  combining  observations 
with  numerical  models  to  improve  forecasting  accuracy.  To  do  this,  observa¬ 
tions  at  a  prescribed  time  are  melded  with  model  predictions  to  estimate  new 
initial  conditions.  These  new  initial  conditions  are  then  in  turn  used  in  the 
next  forecast  run.  The  development  of  assimilation  methods  which  are  used  in 
physical  oceanography  has  always  lagged  behind  meteorology  by  several  years. 
There  is  no  technical  reason  for  this  lag  since  both  fields  use  the  same  general- 
circulation  models  which  have  the  same  level  of  complexity.  A  lack  of  adequate 
oceanographic  data  sets  has  been  cited  as  the  likely  reason  for  the  lag  (Miller 
et  al.,  1997).  The  meteorologists  have  several  orders  of  magnitude  more  data  at 
their  disposal  than  the  oceanographers.  A  sobering  statistic:  the  total  number 
of  recorded  oceanographic  in  situ  observations,  over  the  last  100  years  or  so, 
would  only  fill  2.4  gigabytes. 

Many  approaches  have  been  used  to  assimilate  data  into  comparatively 
complex  ocean  models  using  an  ad  hoc  scheme  in  which  the  observational  data 
“nudges”  the  model  toward  the  measurement  (Malanotte-Rizzoli  and  Holland, 
1986).  The  nudging  methods  fall  into  a  class  similar  to  that  of  direct  data 
insertion.  These  simple  assimilation  methods  do  not  allow  for  the  non-local 
nature  of  the  measurements  themselves  and  neglect  the  off-diagonal  terms  of 
the  model  error  covariance  and  data  covariance.  Gaillard  (1992)  performed  a 
twin  experiment  using  a  two  mode  QG  model.  Point  and  tomographic  data 
were  simulated  and  used  to  reconstruct  the  stream  function  fields  using  sim¬ 
ple  time-independent  inversions.  Although  the  first  baroclinic  mode  was  ade¬ 
quately  resolved,  the  barotropic  mode  was  not. 

1.3  Novel  Aspects  of  This  Work 

This  work  presents  the  results  of  assimilating  tomographic  sum  and  recipro¬ 
cal  travel  time  measurements  into  a  linear  Rossby  Wave  model,  a  non-linear 
Rossby  Wave/Advection  model  and  a  quasi-geostrophic  open  ocean  circulation 
model.  A  Kalman  filter  and  various  sub-optimal  updating  schemes  were  used. 
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Twin  experiments  are  first  used  to  evaluate  the  information  content  in  the 
assimilations.  Actual  AMODE  sum  and  reciprocal  travel  times  (The  AMODE- 
MST  Group,  1994)  are  then  assimilated  using  these  models.  Energy  statistics 
for  the  AMODE  region  are  presented. 

1.4  Thesis  Overview 

The  balance  of  this  thesis  is  divided  into  four  chapters  in  which  we  discuss 
numerical  experiments  involving  point  and  tomographic  data,  and  the  assimi¬ 
lation  of  the  AMODE  travel  time  data  using  three  different  models. 

We  begin  by  discussing  the  theory  of  ocean  acoustic  tomography  and  the 
theory  of  parameter  estimation  using  Kalman  filtering  in  Chapter  2. 

Next,  Chapter  3  outlines  the  models  used  to  fit  the  tomographic  measure¬ 
ments.  These  models  vary  in  their  complexity  and  include  a  linear  Rossby  wave 
model,  with  and  without  advection,  in  addition  to  a  non-linear  open  ocean  cir¬ 
culation  model. 

Chapter  4  discusses  the  results  of  the  parameter  estimation.  It  is  divided 
into  three  major  experiments,  El,  E2,  and  E3.  The  first  portions  of  El  and  E2 
discuss  results  from  numerical  experiments,  and  the  second  portions  of  El  and 
E2  and  E3  discuss  results  of  fitting  real  AMODE  travel  time  data. 

Finally,  in  Chapter  5  we  discuss  the  findings  of  this  thesis  and  offer  oppor¬ 
tunities  for  future  research  directions. 


Figure  1.1:  As  part  of  a  larger  ocean  observational  system,  acoustic  tomogra¬ 
phy  would  be  combined  with  other  forms  of  ocean  measurements,  such  as  ship 
based  CTD  and  satellite  altimeter  data.  Its  integral  and  synoptic  properties 
help  it  to  complement  these  other  measurements. 


Twin  Experiment  Block  Diagram 


Figure  1.2:  Block  diagram  of  a  twin  experiment 


9 


Chapter  2 

THEORY 

This  chapter  discusses  the  theory  of  acoustical  tomography  travel  time  mea¬ 
surements  as  well  as  parameter  estimation.  By  no  means  is  this  an  exhaustive 
discussion;  it  simply  covers  the  main  topics  of  these  subjects  that  are  germane 
to  this  research. 

2.1  Acoustic  Measurements 

As  discussed  in  Chapter  1,  the  sampling  properties  of  tomographic  data  are 
quite  different  from  those  of  point  data.  Ray  theory  provides  a  relatively  simple 
yet  complete  means  of  explaining  acoustic  propagation  for  most  tomographic 
applications.  It  can  generally  be  used  when  the  acoustic  energy  is  not  scattered 
or  absorbed  by  the  ocean’s  surface  or  bottom.  The  AMODE  experiment  was 
conducted  in  deep  enough  water  (~5000  m),  and  the  sources  used  high  enough 
frequency  (250  Hz),  that  ray  theory  can  be  used  to  model  the  acoustic  measure¬ 
ments. 

Using  ray  theory,  we  can  view  the  propagation  of  acoustic  energy  between 
a  source  and  a  receiver  as  being  a  summation  of  multipaths.  A  multipath  is 
an  eigenray  which  satisfies  Snell’s  law  and  intersects  the  source  and  receiver. 
Each  path  is  slightly  different  from  the  next;  an  acoustic  pulse  that  originates 
at  the  source  will  require  a  slightly  different  travel  time  to  complete  the  journey 
to  a  given  receiver. 

The  refraction  of  the  acoustic  paths  is  primarily  controlled  by  the  vertical 
soundspeed  gradient,  and  effects  of  currents  can  be  ignored  when  tracing  the 
rays.  An  average  soundspeed  curve  for  the  AMODE  region  is  shown  in  Figure 
2.1,  Panel  A.  This  soundspeed  curve  was  derived  using  yearly  averaged  Levi- 
tus  hydrographic  data  taken  within  the  AMODE  region  (Levitus,  1982).  The 
soundspeed  minimum  (which  for  this  region  is  at  a  depth  of  about  1200  m)  acts 
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as  a  geometrically  dispersive  wave  guide  or  “sound  channel”.  If  the  source  and 
receiver  pair  are  positioned  near  the  axis  of  this  minimum,  the  ray  paths  will 
“oscillate”  around  the  minimum  axis,  with  turning  points  above  and  below  the 
axis  (See  Panel  B  of  Figure  2.1). 

As  Figure  2.1  shows,  each  ray  has  a  number  of  turning  points.  The  hori¬ 
zontal  distance  between  two  adjacent,  upper  or  lower,  turning  points  is  called 
the  double  loop  length.  Typically,  in  the  AMODE  region  this  distance  is  on  the 
order  of  30  km  for  low  angle  rays  and  about  80  km  for  steeper  angle  rays.  The 
turning  depths  and  the  double  loop  length  of  the  ray  primarily  determine  the 
spatial  sampling  properties  of  the  raypath.  These  sampling  properties  include 
the  depth  resolution  and  limited  horizontal  resolution  of  the  raypath. 

The  travel  time  path  integral  gives  the  time  of  flight  for  an  acoustic  pulse 
along  a  given  ray  path.  This  equation  represents  the  measured  travel  time  T 
of  the  ?;th  raypath  between  the  source  and  receiver  pair. 


Ti  =  f  J-TT T7 — 7T7 — 7 — 7\ - ^  6t  (2.1) 

J  c0(x)  +  dc(x,  t)  +  u(x,  t)  •  r 

where  ds  is  the  incremental  path  length,  t  is  geophysical  time,  c0(x)  is  a  refer¬ 
ence  soundspeed,  8c(x,  t)  is  perturbation  soundspeed  and  u(x,  t)  ■  t  is  the  com¬ 
ponent  of  the  three-dimensional  current  velocity  along  the  raypath.  The  noise 
term  has  been  added  to  represent  instrument  error,  internal  wave  noise, 
etc.  The  travel  times  may  be  summed  and  differenced  between  sources  and 
receivers  in  each  direction  in  order  to  estimate  Sc(x,t )  and  u(x,  t).  After  lin¬ 
earizing  this  equation  about  the  reference  soundspeed  c0(x)  and  forming  sum 
and  differences  of  forward  (+)  and  reciprocal  (-)  perturbation  travel  times  8T, 
we  get: 


6T£  ±  «T 
2 

SJl  -  SIT 

2 


8c(x,  t ) 
c0(x)2 
u(x,  t)  ■ 
Co(x) 


ds  +  e+ 
—ds  +  e_  . 


(2.2) 


An  implicit  assumption  made  in  this  linearization  is  that  the  perturbations  in 
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soundspeed  c0(z)  are  sufficiently  small  that  they  do  not  modify  the  raypath  Tj 
(the  frozen  field  assumption).  For  soundspeed  fields  encountered  in  this  re¬ 
search,  the  non-linear  effects  in  travel  time  are  estimated  to  be  on  the  order 
0.1  ms,  which  is  considerably  less  than  the  noise  level  of  the  data,  e.  The  ad¬ 
ditive  noise  term,  e+,  represents  the  combined  noise  terms  which  may  result 
from  unmodeled  mooring  motion,  clock  drift,  internal  wave  noise,  just  to  name 
a  few.  And  e_  contains  primarily  errors  resulting  only  from  clock  drift. 

In  order  to  invert  the  travel  time  data,  the  time  series  of  measurements 
needs  to  be  carefully  matched  or  aligned  to  the  respective  predicted  eigenrays, 
Tj.  The  modeled  acoustic  rays  were  aligned  with  the  data  by  matching  patterns 
in  the  arrival  times  and  vertical  arrival  angles.  The  end  product  of  this  arrival 
matching  process  identified  about  15  arrivals  for  each  source  and  receiver  pair. 
For  this  study,  three  rays  per  source  and  receiver  pair  were  used,  an  axial, 
steep,  and  an  intermediate  turning  ray.  It  was  determined  that  vertical  resolu¬ 
tion  is  diminished  negligably  by  using  only  a  subset  of  the  ray  paths. 

It  is  convenient  to  think  of  tomography  as  building  an  image,  in  a  sequence, 
of  vertical  and  horizontal  slices  (Kak  and  Slaney,  1988).  The  previous  sections 
discussed  factors  affecting  the  vertical  resolution;  now  we  will  examine  how 
the  map  or  model  may  be  resolved  in  the  horizontal  slice.  The  Fourier  slice 
theorem  tells  us  that  the  Fourier  transform  of  a  projection  is  equivalent  to  the 
Fourier  transform  of  the  image  along  a  given  radial  9. 

Figure  2.2  shows  a  simple  tomographic  geometry  and  diagrams  the  Fourier 
slice  theorem.  It  shows  that  in  order  to  build  up  the  image  in  Fourier  space, 
multiple  projections  at  different  0’s  need  to  be  made.  Once  the  image  is  built  in 
wavenumber  space,  its  equivalent  can  be  obtained  in  physical  space  by  using 
an  inverse  transform. 

In  Chapter  3,  we  will  parameterize  the  soundspeed  perturbation,  5c,  and 
current  variability,  u(x,  t),  as  a  set  of  harmonics  S( k)  (k  is  the  wavenumber 
ordinate,  and  x  the  physical  space  ordinate).  We  will  show  that  the  sum  travel 
time  perturbations  from  Equation  2.2  can  be  re-written  as: 


6T+  = 


f  Zc(z)S( k)  exp[— ik  •  x]dk 


And  the  difference  travel  times  are  given  by: 
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f  Zfaz)S(k)(iVx  exp[-£k  •  x])  •  rdk 
co(z) 


ds  . 


The  complex  amplitudes  S( k)  are  the  Fourier  components  of  the  stream  func¬ 
tion  field. 


2.2  Parameter  Estimation  and  Kalman  Filtering 

The  overall  goal  of  parameter  estimation  is  to  optimally  (what  is  optimal  will 
be  defined  shortly)  fit  a  set  of  model  parameters,  Xk,  to  a  set  of  noisy  mea¬ 
surements.  These  parameters  evolve  between  measurement  times  according  to 
uncertain  dynamical  laws.  Further,  some  or  all  of  the  measurements  and  dy¬ 
namical  transition  functions  may  be  nonlinear,  i.e.,  functions  with  non-constant 
derivatives. 

More  formally,  the  problem  is  to  estimate  the  state  sequence  Xk,  where  xk 
is  a  state  vector  of  dimension  m  for  each  time  step  k.  We  are  given  a  series  of 
vector  valued  measurements  (say  travel  times),  Zk,  as  well  as  an  estimate  of 
the  initial  state,  x0.  Further,  let  n(k)  be  the  dimension  of  the  data  space  for  the 
kth.  time  step. 

We  can  then  relate  the  data  vector,  Zk,  to  the  state  values  using  a  measure¬ 
ment  equation 

Zk  —  hk(xk)  -f~  Vk  ,  (2.3) 

where  hk  represents  a  noiseless  measurement  functional  and  t is  the  noise 
term  which  is  normally  distributed  with  mean  zero  and  variance  Rk  .  We  as¬ 
sume  that  the  functions  h(k)  and  Rk  are  known  a  priori. 

Now  a  transition  equation  can  be  written  which  relates  the  state  estimate 
at  each  time  point  to  the  state  values  at  a  previous  time  point: 

Xk  =  fa{xk-i)  +Wh  ,  (2.4) 

where  fa  represents  the  noiseless  transition  and  Wk  is  normally  distributed 
with  mean  zero  and  variance  Qk .  We  assume  that  the  function  (f>k  and  variance 
Qk  are  known.  The  initial  estimate  of  the  state  x 0  is  related  to  the  true  initial 
by 
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%0  —  £o  +  60  , 

where  e0  is  normally  distributed  with  mean  zero  and  known  variance  P0  . 

We  assume  the  all  the  noise  variances  are  uncorrelated  and  are  positive 
semi-definite  (no  negative  eigenvalues).  The  cost  function,  J,  is  defined  as  the 
negative  log  likelihood  function  for  our  problem  and  is  given  by: 


J(x o,xu...xN)  =  (zo  —  %o)Pq  1(^o  ~  ^o)T  (2.5) 

N 

k= 1 

+  ^(Xk  -  (j)k(xk-i ))Qi1(xk-<l>k(xk-i))T  • 

/c=l 

The  optimal  estimate  for  will  be  the  one  which  minimizes  the  cost  func¬ 
tion.  The  particular  method  used  to  find  the  minimization  is  highly  problem 
dependent.  For  instance,  when  all  of  the  measurement  and  transition  func¬ 
tions  are  affine  (constant  derivative),  an  ordinary  least  squares  method  can 
be  used  which  takes  advantage  of  the  “uncoupled”  nature  of  the  measurements 
and  transition  processes.  This  method  is  known  as  the  Kalman  filter  or  Kalman 
smoother  (the  affine  smoother  is  a  Kalman  filter  with  a  backwards  solution  en¬ 
hancement  portion). 

It  can  be  shown  that  the  Kalman  smoother/filter  in  effect  makes  successive 
transformations  of  the  minimization  problem  for  J  at  each  time  step,  k.  Each 
minimization  corresponds  to  a  step  in  the  filter  or  smoother.  We  will  simply 
state  the  results  of  the  transformation;  the  interested  reader  is  referred  to  Bell 
(1994)  for  the  proof.  The  first  portion  of  the  transformation  is  composed  of  the 
least  squares  update  theorem: 

x  =  x  +  K(z  —  h(x)) 

P  =  (7  —  KH)P 

K  =  PHtS~\  S  =  (HPHt  +  R)  . 
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For  convenience,  we  have  omitted  the  subscripting  for  the  time  step,  and  we 
have  introduced  the  following  terms:  x  is  the  state  estimate  (at  the  kth  time 
step)  derived  from  x  the  forecast  model  estimate.  K  is  the  Kalman  “gain”, 
which  pre-multiplies  the  “residual”,  i.e.,  the  difference  between  z  (the  observa¬ 
tions)  and  h(x)  (the  forecast  estimate  of  the  observations  at  the  k  th  step).  P 
is  the  estimated  covariance  of  the  state  estimate,  I  is  the  identity  matrix,  H  is 
the  gradient  (with  regard  to  x)  of  the  measurement  functional  h,  and  P  is  the 
forecast  covariance  (which  is  discussed  next). 

In  order  to  propagate  the  state  estimate  and  its  covariance  {x  and  P)  to  the 
next  time  step,  we  use  the  transition  step 


P  =  +  Q 

x  =  <t>{x)  , 

where  we  have  introduced  the  new  term  $  which  is  the  mxm  transition  matrix 
composed  of  the  gradient  of  4>  with  regard  to  each  of  the  model  parameters  x. 
This  algorithm  is  collectively  known  as  the  Kalman  filter  and  will  be  used  as 
an  estimation  method  for  the  models  discussed  in  the  following  chapter. 
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Figure  2.1:  Panel  A  shows  the  reference  soundspeed  used  in  the  ray  trace. 
Panel  B  shows  a  sample  of  three  rays  for  a  650  km  path,  and  Panels  C  and  D 
show  the  vertical  eigenmodes  of  soundspeed  and  stream  function  used  in  the 
modeling. 
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Chapter  3 

MODELS 


3.1  Introduction 

This  work  is  designed  after  and  uses  data  collected  from  the  AMODE  experi¬ 
ment,  which  assembled  a  tomographic  array  of  six  acoustic  transceivers  moored 
to  the  ocean’s  bottom  in  a  pentagonal  array  (Dushaw  et  al.,  1996).  The  array 
was  located  between  Puerto  Rico  and  Bermuda  (See  Figure  3.1).  The  AMODE 
experiment  actually  consisted  of  several  experiments.  The  moored  portion  col¬ 
lected  a  long  time  series  of  sum  and  difference  travel  time  data  as  well  as  point 
temperature  data  at  three  depths  along  each  of  the  moorings.  Another  por¬ 
tion  consisted  of  a  high-resolution  moving  ship  tomography  experiment,  which 
collected  one-way  travel  time  data  between  the  moorings  and  a  ship  circum¬ 
navigating  the  perimeter  of  the  array.  A  third  portion  consisted  of  verification 
data  collected  from  CTD  surveys  conducted  in  the  moving  ship  tomography  ex¬ 
periment  and  two  additional  AXBT  surveys.  This  work  will  utilize  a  limited 
portion  of  the  AMODE  data,  that  from  the  moored  transceivers. 

The  ranges  between  the  moorings  were  nominally  350,  410,  and  660  km, 
and  the  resulting  series  of  travel  time  measurements  covered  180-300  days 
with  sampling  occurring  every  4  hours  on  every  4th  day.  All  of  the  travel  time 
data  was  high-pass  filtered  prior  to  analysis  by  B.  Dushaw  in  order  to  remove 
daily  (tidal)  variability. 

The  AMODE  region  lies  to  the  south  of  the  Gulf  Stream  recirculation  zone. 
The  MODE  Group  (1978)  studied  the  area  slightly  to  the  north  and  determined 
that  that  quasi-geostrophic  dynamics  could  be  used  to  model  the  region,  and 
that  geostrophy  and  non-divergence  are  established  to  within  10%.  Figure  3.2 
shows  the  mean  kinetic  energy  (Wyrtki  et  al.,  1976)  for  the  region  as  well  as 
other  primary  oceanographic  features  for  the  area. 
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3.2  Rossby  Wave  Model 

In  this  section,  we  will  develop  the  non-linear  Rossby  wave  model  which  is 
complete  with  advection.  The  linear  version  with  no  advection  follows  trivially 
by  setting  the  advection  to  zero. 

Rossby  waves,  also  called  planetary  waves,  rely  in  part  on  the  earth’s  cur¬ 
vature  (i.e.,  the  variation  of  Coriolis  force  with  north-south  distance)  for  their 
existence.  Rossby  waves  are  solutions  to  the  potential  vorticity  equation  with 
no  forcing.  We  can  separate  the  solution  into  horizontal  and  vertical  parts. 
This  separation  allows  us  to  parameterize  the  solution  for  the  stream  function 
by  using  a  smooth  superposition  of  spectral  functions  in  the  horizontal  and 
quasi-geostrophic  modes,  Zi(z),  in  the  vertical.  These  modes  are  solutions  to 
the  vertical  portion  of  the  quasi-geostrophic  equation  with  flat  bottom  and  free 
surface  boundary  conditions. 

The  horizontal  modes  are  simple  sines  and  cosines,  so  we  can  express  the 
stream  function,  ip(x,z,t),  as  a  summation  of  complex  wave  amplitudes,  %/, 
and  vertical  functions,  Zt(z): 

il>(x,z,t)  =  Y  Zi(z)Sjti  exp(i(k  •  x  -  QJj,it))  ,  (3.1) 

3,1 


where  we  have  summed  over  the  j  horizontal  and  l  vertical  modes.  The  wave 
frequency,  oj,  controls  the  time  dependence  of  the  waves  and  must  satisfy  the 
dispersion  relationship.  The  dispersion  relationship  is  found  by  constraining 
the  waves  to  satisfy  the  vorticity  equation,  which  will  be  discussed  in  Section 
3.3.2.  After  doing  so,  we  find  that  for  a  fluid  with  advection,  ui  has  the  form 


03  = 


k 


-  A uja 


where  K 2  =  k2  + 12,  and  we  have  introduced  A i  which  is  an  eigenvalue  from 
the  solution  to  Poisson’s  equation  shown  in  Equation  3.9.  The  term  A coa  is 
the  portion  of  the  dispersion  relationship  caused  by  the  horizontal  advection 
velocity  for  each  of  the  modes. 

We  will  use  Ui  for  the  advection  velocity  of  the  barotropic  mode,  and  U2  for 
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the  advection  velocity  of  the  baroclinic  mode.  We  find  that  Acoa  for  the  first  two 
modes  is  equal  to: 


Aw0,i 

Awa>2 


-Ui  k 


-2Ui 


U2  •  k eK2 
K2-  \\ 


For  the  linear  Rossby  wave  model,  Ui  and  U2  are  zero  and  we  can  omit  the  A oj 
terms  from  the  dispersion  formula. 

A  diagram  showing  the  wave  parameters,  Sjj,  in  the  complex  plain  is  shown 
in  Figure  3.3.  The  harmonics  Sjj  are  two-dimensional  Fourier  harmonics  in 
a  box,  which  must  contain  the  data  array  but  can  be  any  size.  If  the  domain 
of  the  harmonics  is  the  same  as  the  domain  of  the  data,  then  the  periodicity 
of  the  expansion  functions  will  be  enforced  in  the  inverse.  In  the  experiments 
involving  numerical  simulations  of  data,  the  data  is  periodic  so  that  no  discon¬ 
tinuities  are  present  even  though  the  data  and  model  domain  sizes  coincide.  In 
the  assimilation  runs  using  real  data,  the  model  domain  was  chosen  to  be  two 
times  as  large  as  the  data  domain,  in  order  to  avoid  Gibbs  effects  (Cornuelle 
et  al.,  1989). 

The  a  priori  energy  spectrum  in  the  model,  P,  is  assumed  to  be  “red”,  i.e. 
homogeneous  and  isotropic,  monotonically  decreasing  with  increasing  scalar 
wavenumber  (Wunsch,  1983).  For  these  experiments  we  have  assumed  the 
functional  form  for  the  spectral  power  to  be: 


P(K)  = 


5 


where  K  =  \/k2  + 12,  and  we  use  Ks  as  a  parameter  to  set  the  half  power  point 
“shoulder  point”  for  the  spectral  power.  Figure  3.4  shows  the  spectrum  of  the  a 
priori  covariance. 

The  Kalman  filter  requires  specification  of  the  model  noise  covariance  Q. 
For  this  research,  we  have  chosen  to  use  a  scaled  version  of  the  a  priori  covari¬ 
ance,  P0,  which  relaxes  to  zero  with  a  time  scale,  r  (Howe  et  al.,  1987).  This 
assumes  that  P0  specifies  the  climatological  variances  from  the  reference  state. 
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Q  then  may  be  written  as: 

Q(St)  =  P0(si)  t  <  r 

-  Po  t>T. 

The  time  scale,  r,  was  chosen  to  be  approximately  300  days.  We  will  now  dis¬ 
cuss  the  measurement  model  for  the  travel  time  data. 

3.2.1  Data  model 

The  tomographic  data  which  consists  of  the  sum  and  reciprocal  travel  times 
between  moorings  can  be  expressed  as: 

Zn  (t)  =  '^2Hn,j,iSj,iexv[-iuj,i(U)t)]  (3.2) 

hi 

where  zn(t)  is  the  nth  travel  time  datum  at  the  time  step  t  =  tk.  Hnjyi  is  the 
measurement  functional  that  maps  each  complex  wave  amplitude,  Sjj,  to  a 
particular  travel  time  datum.  For  this  research  we  used  15  source  receiver 
pairs  each  having  3  eigenrays  (45  sum  travel  times  and  45  reciprocal  travel 
times). 

The  gradients  of  Equation  3.2  with  respect  to  the  model  parameters  S  and 

U  =  [iti,  vi,u2,  v2]  are: 


5SU 

(3.3) 

6zn(t) 

Sui 

j,l=  1  jyl—2 

(3.4) 

5zn(t) 
5v  1 

=  e  -mw + E  ~{msuG 

i  1=1  7 /— 9 

(3.5) 

— ikteK 2  c  n 

~  x2  —  A2  ’  x 

j,i= 2  2 

(3.6) 

SU2 

5zn  (t) 
6v2 

— ilttK 2  _ 

—  v 2  •  y  • 

J,*= 2  2 

(3.7) 

Note  that  for  the  linear  Rossby  wave  model,  U  =  0,  and  the  above  equation 
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simplifies  to: 

SZn(t)  _  rr 

SShl  n'hl  ■ 

For  the  non-linear  case,  the  additional  terms  in  Equation  3.3  must  be  evalu¬ 
ated.  We  will  now  discuss  the  non-linear  circulation  model,  after  which  we  will 
return  to  discussion  of  the  Rossby  wave  models  again.  The  assimilation  results 
are  given  in  Chapter  4. 

3.3  Non-Linear  Open-Ocean  Model 

3.3.1  Introduction 

Within  the  last  few  years,  the  ocean  modeling  community  has  made  great 
strides  and  advances.  Their  efforts  have  led  to  an  unprecedented  variety  of 
complexity  in  ocean  circulation  models  (Holland  and  Capotondi,  1996).  Some 
models,  such  as  those  that  couple  the  circulation  of  the  atmosphere  and  global 
oceans,  are  very  complex  and  require  a  great  deal  of  computational  resources. 
Others,  such  as  single-layer  gravity-wave  models,  are  more  simplistic  and  re¬ 
quire  relatively  meager  computational  resources.  Quasi-geostrophic  (QG)  mod¬ 
els  fall  at  about  the  midpoint  of  this  spectrum.  QG  models  use  a  simplified  set  of 
equations  to  describe  the  motion  of  large-scale  and  meso-scale  features.  These 
motions  result  from  slight  force  imbalances  that  exist  between  the  geostrophic 
pressure  gradients  and  coriolis  forces  (Apel,  1987).  The  pressure  gradients  are 
directly  related  to  horizontal  gradients  in  the  water  density. 

This  study  uses  a  QG  model  known  as  the  Harvard  Open  Ocean  Model  or 
HOOM  (Miller  et  al.,  1983).  HOOM  belongs  to  a  class  of  circulation  models 
known  as  regional  or  “open  ocean”  circulation  models.  Regional  models  do  not 
use  explicit  boundaries  or  coast  lines  to  contain  the  circulation  model.  Implicit 
or  computational  boundaries  are  used  to  communicate  or  interact  with  the  flow 
regions  outside  the  model  domain.  Since  the  modeling  domain  need  not  encom¬ 
pass  an  entire  ocean  basin,  the  regional  model’s  horizontal  resolution  may  be 
fine  enough  to  resolve  mesoscale  features  (order  50  km).  This  boundary  com¬ 
munication  involves  specifying  boundary  conditions  along  the  perimeter  of  the 
model.  The  boundary  information  is  only  specified  along  inflow  segments  of  the 
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boundary;  the  outflow  is  prescribed  by  the  interior  flow  field  of  the  model. 

Many  researchers  have  used  HOOM  in  a  variety  of  assimilation  and  mod¬ 
eling  studies.  Rienecker  and  Miller  (1991)  used  the  model  to  investigate  the 
California  Current,  while  Pinardi  and  Robinson  (1987)  used  the  model  in  the 
POLYMODE  region.  Haidvogel  et  al.  (1980)  originally  developed  HOOM  as  a 
barotropic  model,  and  HOOM  was  later  configured  into  a  baroclinic  version  by 
Miller  et  al.  (1983). 


3.3.2  Theory 

The  HOOM  model  integrates  the  quasi-geostrophic  equations  of  motion  on  a 
beta  plane.  The  potential  vorticity,  (,  with  forcing  F  is  integrated  level  by  level 
according  to  the  equation 


d_ 

dt 


+  60J(lp,) 


<  +  f,t  =  F- 


(3.8) 


The  stream  function,  ip,  is  solved  as  a  function  of  the  potential  vorticity  using 
the  Poisson  equation: 

V2ip  +  T2(oipz)z  =  C  ,  (3.9) 

where  J  is  the  Jacobian  operator  defined  as 


J(A,  B)  = 


6A5B  _  5B_5A 
6x  8y  5x  5y  ’ 


(3.10) 
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and  the  following  terms  are  defined  as 


e/5  = 

P  = 
r  = 

a(z)  = 
N{z )2  = 


A 

/3d2’ 

20 

—  cos  0, 
a 

foD 


N0hT ’ 
iV02 


iV(2)2’ 
9_Sp 
Po  Sz' 


(3.11) 

(3.12) 

(3.13) 

(3.14) 


Table  3.1  describes  all  the  variables  introduced  in  the  above  equations,  and 
lists  the  numerical  values  for  the  model  parameters. 

The  stream  function,  ip,  is  directly  related  to  the  pressure  anomaly  using 
the  relationship 

P  =  PoM  •  (3.15) 

The  density  anomaly,  S,  is  often  used  as  a  variable  that  is  assimilated.  It  may 
be  derived  from  ip  by  applying  the  hydrostatic  equation.  We  merely  state  the 
result  as 

5  =  r2o-^  .  (3.16) 

dz 

The  numerical  differencing  of  the  vorticity  equation,  Equation  3.8,  uses  an 
Adams-Bashforth  method  in  time  and  finite  element  method  in  space.  These 
methods  are  respectively  second  and  fourth-order  accurate.  Most  of  the  numer¬ 
ical  details  of  HOOM  are  discussed  by  Haidvogel  et  al.  (1980).  We  present  a 
schematic  layout  of  the  model  in  Figure  3.5. 


3.3.3  Computational  parameters 

Our  model  configuration  consists  of  six  levels,  each  placed  at  depths  corre¬ 
sponding  to  the  zero  crossings  of  the  6th  baroclinic  mode.  In  order  to  totally 
encompass  the  tomographic  array,  and  minimize  harmonic  boundary  effects, 
we  chose  to  use  a  1280  km  square  (97x97  grid)  for  the  control  ocean.  And 
we  chose  the  assimilation  domain  to  be  825  km  (61x61),  so  that  the  boundary 
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would  be  updated  by  the  tomographic  array. 

Table  3.1  lists  the  numerical  parameters  used  in  the  QG  model.  The  best 
values  for  the  AMODE  region  were  used.  Miller  et  al.  (1983)  discusses  the 
vertical  structure  of  the  baroclinic  version  of  the  model,  and  describes  the  so¬ 
lution  methods  for  Poisson’s  equation  3.9.  Figure  3.5  shows  the  ip  modes  and 
the  depths  that  correspond  the  model’s  vertical  levels.  Figure  3.6  shows  the 
bathymetric  data  used  in  the  model  simulations. 

The  buoyancy  profile  used  to  construct  the  model’s  vertical  stratification  is 
shown  along  with  other  hydrographic  data  in  Figure  3.7.  These  profiles  are 
spatially  averaged  profiles,  estimated  from  annual  Levitus  data. 

The  integration  of  Equation  3.8  requires  initial  and  boundary  conditions. 
For  these  boundary  conditions,  we  use  the  Charney-Fjortoft-von  Neuman 
(CFvN)  boundary  conditions  (Charney  et  al.,  1950).  The  CFvN  conditions  dic¬ 
tate  that  we  specify  the  normal  velocity  everywhere  on  the  boundary  and  spec¬ 
ify  the  potential  vorticity  where  the  normal  velocity  is  directed  inward. 

Miller  and  Bennett  (1988)  demonstrated  some  of  the  difficulties  that  may 
arise  near  boundary  points  where  the  imposed  normal  velocity  changes  sign 
(the  tangential  flow  points).  Near  these  points,  inconsistencies  develop.  Ul¬ 
timately,  these  ill-posed  boundary  inconsistencies  advect  into  the  interior  of 
the  model  and  contaminate  the  flow  field.  This  ill-posed  nature  is  significant 
because  the  vorticity  gradients  grow  in  time  and  without  bound.  The  incom¬ 
patible  boundary  conditions  are  unavoidable  in  practice  and  will  cause  the 
model  to  fail  on  time  scales  comparable  to  those  that  characterize  the  boundary 
conditions.  These  problems  highlight  the  definite  need  for  data  assimilation 
in  practical  ocean  forecasting.  New  information  must  be  introduced  into  the 
model’s  interior  at  intervals  that  are  short  compared  to  the  forcing  time  scales. 
Whereas  previous  studies  used  four  discrete  waves  to  construct  regional  bound¬ 
ary  conditions  (see  McWilliams  and  Flierl  (1976)),  this  research  uses  a  complete 
spectrum  of  waves. 

Bottom  bathymetry  plays  a  paramount  role  in  realistic  simulation  of  the 
ocean  flow  field,  primarily  because  of  its  influence  on  vortex  stretching.  The 
bottom  bathymetry  for  the  model  is  parameterized  in  terms  of  a  nominal  depth 
(H0)  and  a  bottom  perturbation,  h(x,  y),  that  is  scaled  by  the  Rossby  number,  e 
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(not  to  be  confused  with  data  error).  For  the  bottom  boundary,  the  bathymetry, 
z,  is  given  by 


z  =  —Ho  +  eh(x,y)  .  (3.17) 

Equation  3.17  expresses  the  topographic  variation  in  terms  of  h,  but  only  in¬ 
cludes  topographic  variations  that  are  consistent  with  the  Rossby  wave  expan¬ 
sion.  The  bottom  topography  interacts  with  the  flow  field  by  inducing  a  vertical 
velocity,  w: 

D  ~ 

w  =  — cn/>2  =  uVh  =  J(ip,  h)  (3.18) 

6 

Time-stepping  the  bottom  velocity  requires  the  total  derivative  of  Equation 
3.18  to  be  re-written  in  terms  of  a  partial  derivative: 

Sd  /s 

-  (aipz)  =  h)  -  uV(cr^)  =  Jty,  h  -  a^Zb)  (3.19) 

The  last  term  in  the  Jacobian  operation  of  Equation  3.19,  aipZb,  is  typically 
called  “bottom  density”.  Equation  3.19  shows  that  the  time  evolution  of  the 
bottom  density  is  made  up  of  both  a  “topographic  term”  and  an  advection  term. 
This  advection  term  keeps  track  of  bottom  induced  effects  from  previous  time 
steps.  Equation  3.19  has  a  form  identical  to  that  of  the  vorticity  equation, 
3.8.  Numerically,  the  model  solves  for  both  the  time  evolution  of  vorticity  and 
bottom  density  in  the  same  way.  Like  the  vorticity  solution,  the  solution  for  the 
bottom  density  requires  boundary  conditions  as  well  (they  must  be  specified 
consistent  with  the  CFvN  conditions). 

3.4  Model  Interpolation  and  Soundspeed 

Tomographic  travel  times  measure  integral  properties  of  soundspeed  anomaly 
(Sc)  and  current  velocity,  u,  according  to  Equation  2.2.  This  poses  two  problems 
when  trying  to  use  these  travel  times  with  numerical  circulation  models: 

1.  The  tomography  data  does  not  directly  measure  properties  of  the  circu¬ 
lation  model  (ip,  Of  so  a  transformation  procedure  must  be  used  to  back 
these  parameters  out  and  also  interpolate  the  model  between  the  layer 
depths. 
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2.  HOOM,  like  most  other  circulation  models,  discretely  represents  the  ver¬ 
tical  structure  of  the  stream  function  and  vorticity  fields  at  a  limited  num¬ 
ber  of  depths.  Because  the  model  specifies  the  stream  function  and  vor¬ 
ticity  only  at  a  few  discrete  levels,  interpolation  methods  must  be  used  to 
get  a  smooth  continuous  soundspeed  field.  Modes  or  empirical  orthogonal 
functions  (which  are  fitted  with  smoothing  splines)  can  be  used  to  con¬ 
struct  the  interpolation  functions.  This  method  has  been  used  by  many 
others,  but  we  will  give  a  succinct  review  of  modes;  see,  for  example,  Cor- 
nuelle  (1983)  as  well  as  Gaillard  (1992). 

The  goal  of  this  section  is  to  demonstrate  a  method  for  interpolating  F(z)  that 
is  based  on  a  normal  mode  expansion  for  the  stream  function  ip.  The  first  in¬ 
terpolation  method  that  we  will  discuss  is  based  upon  the  so-called  “dynamic 
modes”.  The  second  method  uses  a  statistical  basis  known  as  EOFs  (empiri¬ 
cal  orthogonal  functions).  We  will  first  discuss  the  formulation  of  the  dynamic 
modes. 

We  begin  by  assuming  that  the  soundspeed  field,  c(x,  t),  and  density  field, 
p(x,  t),  can  be  expressed  as  first  order  perturbation  fields: 

c(x,t)  =  5c(x,t)  +  co(z)  (3.20) 

p(x,i)  =  Sp(x.,t)  +  p0(z) ,  (3.21) 

The  perturbation  fields,  5c,  5p,  can  be  further  separated,  not  out  of  necessity, 
rather  for  convenience,  into  a  vertical  portion  and  a  horizontal  portion  of  the 
form  5c,8p(x.,t )  =  Fh(x,  y)Fz{z).  The  vertical  modes,  Fz(z),  are  based  on  solu¬ 
tions  to  Poisson’s  equation  3.9  which  expresses  the  relationship  between  the 
stream  function  variable  ip  and  vorticity  £•  Poisson’s  equation,  Equation  3.9,  is 
separable,  that  is  f/>(x)  =  F(x,  y)Z(z).  The  stream  function  may  be  written 

N 

ip(x,t)  =  Y^ki(x,y,t)Z(z) ,  (3.22) 

i= 1 

where  N  is  less  than  or  equal  to  the  number  of  model  layers.  The  expansion 
for  the  vorticity,  £,  uses  the  same  vertical  basis  set  as  does  ip.  This  follows  from 
the  definition  £  =  W2ip. 
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After  Equation  3.22  is  used  in  Poisson’s  equation  (3.9),  the  vertical  portion 
of  the  stream  function  equation  may  be  written 

A2  7 

r2a-f  =  -XZ,  (3.23) 

52z 

with  boundary  conditions  specified  at  the  surface  z  —  0  (although  none  are 
specified  in  this  research)  and  at  the  nominal  model  bottom  z  =  —H  given  by 

^(0)  =  ¥(-H)  =  0  .  (3.24) 

dz  bz 


We  numerically  solve  Equation  3.23  with  these  boundary  conditions  using  an 
integrator  and  shooting  method  in  order  to  obtain  the  modes  Z(z). 

The  soundspeed,  5c,  can  be  obtained  using  the  gradient  of  the  potential 
sound  speed,  Cp,  using  the  following  relationship: 


fo  Sip  5Cp 
N2(z)  5z  5z 


(3.25) 


Note  that  the  gradient  of  the  potential  soundspeed  is  not  specified  by  the  cir¬ 
culation  model  and  must  be  obtained  from  regional  hydrology  (see  Figure  3.7). 
After  accounting  for  various  scaling  constants  used  in  the  model  and  listed 
in  Table  3.1,  we  can  finally  write  the  expression  for  the  soundspeed  modes  in 
terms  of  the  stream  function  modes  as 


Z{z) 


Vo H\  2  6^50, 
f0d  )  °  5z  5z 


3.4.1  Comparisons  Between  Statistical  and  Dynamic  Modes 

The  statistical  modes  (typically  called  empirical  orthogonal  functions,  EOFs) 
are  based  solely  on  statistics  of  the  data  alone.  No  physical  model  or  other  con¬ 
straints  are  used.  They  allow  a  more  efficient  representation  of  the  field  than 
do  dynamic  modes.  However,  they  require  some  type  of  spline  interpolation 
between  the  model’s  layers  which  is  required  in  the  raypath  estimation. 

In  order  to  evaluate  the  trade-offs  between  using  the  dynamic  modes  versus 
the  statistical  modes  the  regional  climatology  for  the  AMODE  region  was  simu¬ 
lated  using  the  QG  model.  In  this  simulation,  a  ten  year  climatological  control 
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run  was  made  using  HOOM  on  a  97x97x6  model  domain.  The  parameters  and 
bottom  topography  are  shown  in  Table  3.1  and  Figure  3.6.  These  parameters 
are  taken  as  being  indicative  of  the  AMODE  region. 

We  then  sampled  this  climatological  simulation  and  estimated  the  amount 
of  variance  that  each  of  the  basis  sets  (dynamical  and  EOF)  explained.  Table 
3.2  shows  the  percentage  of  the  model  variance  explained  by  the  dynamic  and 
EOF  modes.  The  first  two  dynamic/EOF  modes  accounted  for  98.4%/99.2%  of 
the  overall  variance.  Tables  3.3  and  3.4  show  these  variances  for  each  level  in 
the  QG  model  that  is  explained  using  the  first  two  modes.  As  can  be  seen,  both 
sets  of  modes  explain  the  variance  in  the  second  (400  m)  level  the  best  (because 
the  stream  function  and  soundspeed  at  this  depth  are  fairly  uniform).  Finally, 
Table  3.5  shows  the  coefficients  of  correlation  between  the  two  basis  sets.  As 
could  be  expected,  each  of  the  first  EOF  modes  most  closely  correlates  with 
the  first  dynamic  mode,  etc.  Based  on  these  results,  we  believe  that  the  errors 
incurred  from  using  a  truncated  set  of  two  dynamical  modes  to  reconstruct  the 
QG  model  are  negligible.  And  it  was  not  necessary  to  use  EOF  modes  for  the 
purposes  of  this  study. 
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Table  3.1:  Computational  model  parameters 


Symbol 

Name 

Value 

U0 

velocity  scale 

6.01  m  s-1 

to 

time  scale 

1.77590106 sec 

lo 

latitude 

25° 

fo 

Coriolis 

6.1467310-5 

d 

x  scale 

27.2155  km 

H 

vertical  scale 

600.0  m 

L 

vertical  scale 

27.558  (no  units) 

a 

f3  Rossby  number 

00 

t-H 

11 

o 

St 

time  step 

0.00395  days 

r2 

strat.  scale 

f$d2/N$H2  =  0.383912 

N't 

mid-thermocline  buoyancy 

2.02510-5  sec"2 

N2 

climatological  buoyancy 

a(z) 

stability  profile 

N2/N2(z) 

</> 

stream  function  or  pressure 

c 

dynamic  vorticity 

R 

relative  vorticity 

v2v» 

T 

thermal  vorticity 

F2(o"lpz)z 

F 

Shapiro  filter 

K 

kinetic  energy 

u2  +  v2/2 

A 

avail,  pot.  energy 

r2aip2/2 

Table  3.2:  Percentage  of  the  stream  function  variance  explained  by  each  dy¬ 
namic  and  empirical  mode. 


Mode  # 

Dynamic 

EOF 

1 

85.8 

83.7 

2 

12.6 

15.5 

3 

1.36 

0.07 

4 

0.04 

0.02 

5 

0.10 

0.01 

6 

0.10 

<0.01 
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Figure  3.1:  Regional  map  showing  the  AMODE  region.  The  AMODE  array  is 
the  pentagonal  array;  acoustic  transceivers  were  on  each  mooring.  The  QG 
circulation  model  domain  is  indicated  with  the  larger  dashed  box. 


Table  3.3:  Fraction  of  the  variance  of  the  stream  function  fields  which  is  de¬ 
scribed  by  each  dynamical  mode  according  to  depth. 


Depth  level 

DM1 

DM2 

Residuals 

1 

59.4 

37.9 

2.8 

2 

71.3 

28.0 

0.7 

3 

2.4 

4 

ID 

2.5 

5 

94.2 

5.6 

2.6 

6 

87.3 

DSD 

1.3 
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Figure  3.2:  Surface  eddy  kinetic  energy  of  the  mean  flow  for  the  North  Atlantic 
Region,  from  Wyrtkiet  al.  (1976). 


Figure  3.3:  Complex  wave  parameters  Sjj  which  are  used  in  the  Rossby  wave 
model.  Two  modes,  a  barotropic  and  first  baroclinic,  are  used.  Each  mode  has  a 
set  of  40  complex  wave  amplitudes  plus  the  mean  level,  for  a  total  or  162  model 
parameters  (166  with  the  advection  parameters). 
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Spectral  Power 


Figure  3.4:  Log  plot  of  the  a  priori  model  spectrum. 


Table  3.4:  Fraction  of  the  variance  of  the  stream  function  fields  which  is  de¬ 
scribed  by  each  empirical  mode  according  to  depth. 


— 

iDtoiii 

1 

81.3 

1.5 

2 

87.1 

3 

95.6 

2.4 

4 

87.6 

11.4 

5 

80.0 

19.9 

INMOSH 

6 

72.4 

25.6 

| 

33 


Lateral  Boundary 
Conditions 


Z=1100m 


Surface  Forcing  (not  used) 


Level  1  (Stream  Function,  Vorticity) 
Z=100m 


Z=2150m 


Z=3720m 
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Figure  3.5:  Schematic  of  the  open  ocean  QG  circulation  model. 


Table  3.5:  Cross  correlation  energies  between  the  empirical  and  the  dynamical 
modes. 


DM1 

DM2 

DM3 

DM4 

DM5 

DM6 

EMI 

0.701 

0.015 

<0.001 

<0.001 

<0.001 

<0.001 

0.269 

0.773 

0.050 

0.002 

0.004 

<0.001 

0.180 

0.851 

<0.001 

0.034 

0.006 

1MM1 

0.010 

0.035 

0.947 

0.059 

0.012 

EM  5 

0.002 

0.007 

0.004 

0.569 

0.323 

0.003 

0.057 

0.046 

0.334 

0.659 
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Figure  3.6:  Part  A  shows  the  bathymetry  based  on  the  ET0P05  database 
and  resampled  at  the  model  grid  interval.  Part  B  shows  a  contour  plot  of  the 
bathymetry  data  used  in  the  QG  model  simulation.  The  data  was  taken  from  a 
subset  of  the  ET0P05  data  set  and  filtered  harmonically  to  be  consistent  with 
the  QG  dynamics. 
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Figure  3.7:  Hydrographic  data  used  to  construct  the  model.  Plots  of  tempera¬ 
ture,  salinity,  and  soundspeed  modes  as  a  function  of  depth. 
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Chapter  4 

RESULTS 

In  this  chapter,  we  present  results  from  three  different  experiments,  El, 
E2,  and  E3.  In  the  first  two  experiments,  we  will  first  run  a  numerical  twin 
experiment  in  order  to  validate  the  assimilation  method  and  evaluate  how  well 
the  measurements  constrain  the  model. 

In  part,  the  evaluation  is  based  on  how  well  the  measurements  are  able 
to  estimate  the  spectrum  of  the  barotropic  and  first  baroclinic  planetary  wave 
field.  Also,  we  will  assess  the  root  mean  square  differences  of  the  fields,  aver¬ 
aged  over  the  entire  model  domain.  The  goodness  of  the  fits  is  evaluated  using 
three  criteria:  1)  the  error  weighted  sum  of  the  final  residuals  (an  approximate 
X2  distribution),  2)  the  correlation  coefficient,  C,  between  the  observed  and  pre¬ 
dicted  signal  energy,  and  3)  the  percentage  of  the  signal  variance  resolved  by 
the  model,  E.  As  used  here,  the  x2  distribution  should  have  mean  one  with 
variance  equal  to  y^2  /N  (where  N  is  the  number  of  observations  used  in  the 
fit). 

In  El  we  conduct  a  twin  experiment  to  show  the  effectiveness  of  assimilat¬ 
ing  point  and  tomographic  data  using  the  linear  Rossby  wave  model.  After  pre¬ 
senting  the  results  of  these  twin  experiments,  we  assimilate  the  actual  travel 
time  data. 

Next,  in  E2,  mean  large  scale  advection  is  added  to  the  Rossby  wave  model. 
Again,  a  twin  experiment  is  run,  this  time  using  simulated  tomographic  data. 
AMODE  data  is  then  assimilated  and  used  to  estimate  the  wave  field  and  the 
large  scale  mean  flow. 

Finally,  in  experiment  E3,  we  use  a  bootstrap  approach  to  assimilate  the 
AMODE  data  into  the  non-linear  QG  model.  We  call  it  bootstrap  because  we 
will  use  the  Rossby  wave  spectrum  estimated  in  El  to  force  the  boundary  of 
the  HOOM  model  and  re-use  the  AMODE  travel  times  to  update  the  interior 
of  the  model.  The  assimilation  method  will  be  similar  to  that  used  in  El  and 
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E2,  but  no  transition  step  is  applied  to  the  forecast  covariance  update.  This 
form  of  updating  is  called  01  (Optimal  Interpolation),  and  is  numerically  much 
more  efficient  than  the  full  Kalman  update.  It  omits  the  covariance  transition 
between  assimilation  time  steps,  which  contains  order  M2  operations,  where  M 
is  the  number  of  model  parameters. 

4.1  Experiment  El  -  Linear  Rossby  Waves 

In  this  experiment,  we  fit  the  sum  and  difference  travel  times  as  well  as  point 
data  to  a  linear  Rossby  wave  model  which  consists  of  a  barotropic  mode  and 
the  first  baroclinic  mode.  A  twin  experiment  is  run,  first  using  simulated  data, 
then  using  simulated  noisy  data.  After  the  results  of  the  twin  experiments  are 
evaluated,  the  actual  AMODE  data  is  used  to  estimate  the  wave  field. 

In  the  twin  experiments,  simulated  data  are  used  to  update  the  model  at 
1-day  increments.  We  carry  the  assimilations  out  for  a  period  of  300  days, 
which  should  allow  all  components  of  the  wave  field  to  be  adequately  sampled 
in  time  and  space. 

At  each  assimilation  time  step,  the  model  is  updated  using  three  sum  and 
difference  travel  times.  The  model  noise,  Q,  is  prescribed  as  discussed  in  Sec¬ 
tion  3.2.  The  results  of  the  tomographic  data  assimilation  are  shown  in  Figure 
4.1.  We  added  random  noise  to  the  sum  and  reciprocal  travel  times.  This  noise 
corresponded  to  that  which  was  typical  for  the  AMODE  data  and  was  reduced 
by  \/5  in  order  to  partially  account  for  using  3  rays  per  source  receiver  rather 
than  15  which  where  typically  observed.  This  equivalent  error  level  for  the  sum 
and  reciprocal  travel  times  than  amounted  to  ±4.5  ms  and  ±1.3  ms  respectively. 
Panel  A  shows  the  RMS  error  between  the  control  ocean  and  the  assimilation 
estimated  at  each  time  step  (after  update)  in  cyan.  As  can  be  seen,  it  slowly 
decreases  monotonically  during  the  run.  Panel  A  also  shows  the  RMS  error 
and  normalized  x2  as  a  function  of  time,  x2  was  normalized  by  the  number 
of  data  used  in  the  update.  In  this  simulation  it  is  initially  near  one,  but  de¬ 
creases  slightly  to  about  0.7  near  the  end  of  the  300  day  time  period.  This  slight 
decrease  is  due  to  Q  being  specified  slightly  too  large. 

Also  displayed  in  Panel  A  are  the  RMS  variance  in  the  predicted  data  (black) 
as  well  as  the  RMS  of  the  residuals  after/before  the  update  at  each  time  step 
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(blue/yellow).  It  is  useful  in  evaluating  the  success  of  the  assimilation  process 
to  compare  the  data  RMS  level  with  that  of  the  residuals.  The  better  the  fit, 
the  greater  the  reduction  in  variance  will  be. 

Panel  B  shows  the  final  assimilation  map  of  the  baroclinic  mode  at  the  end 
of  the  300  day  run.  The  maps  of  the  eddy  fields  are  virtually  identical,  with 
slight  discrepancies  near  the  edges  of  the  model  122  x  122  domain,  far  from  the 
tomographic  array. 

Panels  C  and  D  show  of  the  error  spectra  as  a  function  of  time  for  the 
barotropic  and  baroclinic  modes  respectively.  Keep  in  mind  that  these  are  error 
estimates  for  a  single  forecasting  run.  An  ensemble  average  of  numerous  runs 
would  be  necessary  to  statistically  estimate  a  more  rigorous  error  field.  But, 
this  single  estimate  is  useful  for  this  discussion. 

For  comparison,  Panels  E  and  D  present  the  analysis  error  estimates  from 
the  assimilation  given  by  the  diagonal  elements  of  P.  As  can  be  seen,  certain 
wavenumbers  are  not  well  resolved  by  the  array.  And  the  errors  in  these  par¬ 
ticular  wavenumber  estimates  persist  in  time,  so  that  the  estimate  at  the  end 
of  the  assimilation  for  these  waves  is  little  better  than  the  initial  estimate. 

Next,  another  twin  experiment  was  performed  using  the  linear  Rossby  wave 
model  with  point  data.  In  this  experiment,  we  will  only  estimate  the  first  baro¬ 
clinic  mode  and  we  will  extent  the  assimilation  period  out  to  600  days  (in  order 
to  sample  the  wave  field  several  times  over).  Also,  we  will  extend  the  wavenum¬ 
ber  domain  out  to  8  cycles  per  megameter  from  the  3.1  cycles  per  megameter 
used  previously. 

The  point  data  consists  of  an  array  of  16  stations  evenly  spaced  throughout 
the  model  domain.  Each  station  sampled  the  water  column  at  three  depths 
in  order  sample  the  first  baroclinic  mode  the  same  number  of  times  as  the 
tomographic  data.  These  point  measurements  could  be  thought  of  as  being 
equivalent  to  zero  length  tomographic  array  sampling  point  travel  time  or  tem¬ 
perature  measurements.  In  order  to  make  the  point  and  tomographic  measure¬ 
ments  as  similar  as  possible,  we  will  only  use  3  sum  travel  times  between  each 
source  and  receiver. 

A  point  data  error  equivalent  to  90  ms  was  assigned  to  the  point  data  in 
order  for  it  to  have  an  equal  weighting  with  the  tomographic  measurements. 
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This  error  level  resulted  in  both  the  tomographic  and  point  data  reducing  the 
error  variance  at  the  end  of  the  600  day  assimilation  period  to  the  same  average 
or  RMS  level  (over  all  wavenumbers).  We  will  compare  the  two  measurements 
by  plotting  the  error  variance  in  the  wavenumber  estimates  for  both  data  sets 
at  a  number  of  times  during  the  assimilation  run  (see  Figure  4.2).  The  point 
data’s  error  variance  is  shown  with  dashed  lines,  and  the  tomographic  data’s 
with  solid  lines.  As  more  data  is  assimilated,  the  point  array’s  estimate  of  all 
wavenumbers  decreases  to  its  final  level  (marked  with  day  600),  while  for  the 
tomographic  data  it  decreases  for  the  low  wavenumbers  only.  Although  both 
data  sets  reduce  the  average  error  in  the  wavenumber  estimates  to  the  same 
level,  the  tomographic  data’s  final  estimate  for  the  smaller  wavenumbers  is 
better  than  that  of  the  point  data. 

Now,  we  will  discuss  results  of  fitting  the  actual  AMODE  data.  Figure  4.3, 
shows  the  wave  parameters  estimated  using  the  real  AMODE  data.  Like  the 
simulations,  only  six  travel  times  (3  sum/3  difference)  per  source/receiver  pair 
were  used  in  the  fit.  In  this  case  the  measurement  error  was  not  deceased  by 
the  \/5  factor;  instead  the  actual  measurement  error  of  the  data  itself  was  used. 

Panel  A,  as  in  the  previous  experiment,  displays  measures  of  how  well  the 
model  is  fitting  the  AMODE  data.  The  x2  estimates  of  the  fit  are  reason¬ 
able,  as  well  as  the  variance  of  the  residuals  and  the  innovations  (differences 
between  forecast  data  and  measured  data).  The  high  correlation  coefficient, 
C  =  0.84,  and  fairly  large  percentage  of  the  data  variance  explained  by  the 
model,  E  =  60.7%,  indicate  that  the  linear  Rossby  wave  model  fits  the  data  rea¬ 
sonably  well.  Panel  B  shows  a  map  of  the  final  baroclinic  stream  function  field. 
Panel  C  shows  that  the  barotropic  mode  has  some  energy  at  higher  wavenum¬ 
bers.  The  baroclinic  mode  shown  in  Panel  D  has  predominantly  low  wavenum¬ 
ber  energies. 

Although  the  fit  to  the  linear  model  is  quite  reasonable,  the  model  does  not 
admit  any  large  scale  advection.  The  next  section  discusses  results  of  fitting 
the  Rossby  wave  field  with  advection. 
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4.2  Experiment  E2  -  Rossby  Waves  with  Advection 

In  Experiment  E2  we  essentially  repeat  the  previous  experiments  with  the  ad¬ 
dition  of  advection  to  the  model.  The  advection  velocity  causes  Doppler  shifts 
in  the  frequencies  of  the  Rossby  waves.  Our  model  now  consists  of  the  ampli¬ 
tude  and  phases  of  the  Rossby  waves,  as  well  as  the  components  of  the  advec- 
tive  velocity  U.  Figure  4.4  shows  the  twin  experiment  results  using  simulated 
tomographic  data  with  the  travel  time  noise  (reduced  by  \/5  as  discussed  in 
Section  4.1)  added  to  the  sum/difference  travel  time  measurements. 

Again,  Panel  A  displays  statistics  of  the  data  fit.  The  x2  values  from  the  fit 
and  the  variance  reductions  are  about  the  same  as  those  of  the  linear  Rossby 
waves  of  experiment  El.  The  correlation  coefficient  and  the  error  variance 
explanation  are  high  for  the  model.  The  true  advection  velocities  for  the  twin 
experiment  are  Uf,c  =  —0.02  m/s,  Ubt  =  —0.02,  m/s. 

The  light  blue  line  in  Panel  A  shows  the  RMS  error  in  the  estimated 
barotropic  advection  speed.  It  does  not  decrease  as  more  travel  times  are  as¬ 
similated,  indicating  that  the  data  is  not  able  to  resolve  the  barotropic  advec¬ 
tion  velocity.  However,  as  shown  by  the  light  yellow  line,  the  estimate  of  the 
barotropic  advection  velocity  does  converge  to  the  true  value.  Panel  B  shows 
the  final  estimated  baroclinic  stream  function  map.  It  is  very  similar  to  that  of 
El.  The  errors  in  the  spectral  estimates  of  the  stream  function  shown  in  Panels 
C  and  D  are  also  similar  to  those  found  in  El. 

The  results  of  fitting  the  AMODE  data  are  shown  in  Figure  4.5.  Panel  A 
show  that  the  x2  estimates  of  the  fit  are  reasonably  close  to  unity  as  in  Experi¬ 
ment  El,  and  that  the  variance  reductions  are  fairly  good. 

Panels  C  and  D  also  show  that  the  estimated  spectrum  of  Rossby  waves  is 
nearly  identical  to  that  of  the  non-advective  model.  Panel  B  shows  the  esti¬ 
mated  advection  velocities  for  both  modes.  The  estimates  for  the  advection  of 
both  modes  are  not  statistically  different  from  zero. 

The  correlation  of  the  fit  and  the  variance  reduction  (0.85/62.15%)  are  virtu¬ 
ally  identical  to  the  linear  Rossby  wave  fit.  Given  that  the  twin  simulations  in¬ 
dicated  that  the  tomographic  data  were  able  to  resolve  the  baroclinic  advection 
rates,  it  is  reasonable  to  suggest  that  baroclinic  advection  is  not  a  significant 
factor  in  the  AMODE  region. 
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4.3  Experiment  E3  -  Non-linear  HOOM  model 

In  Experiment  E3,  the  travel  time  data  were  assimilated  into  the  non-linear 
QG  model.  The  necessary  boundary  conditions  for  the  circulation  model  were 
supplied  using  the  Rossby  wave  field  estimated  with  the  Kalman  filter  in  Ex¬ 
periment  El.  The  boundary  conditions  were  applied  to  all  six  levels  of  the 
model  at  an  update  rate  of  once  per  day. 

Before  assimilating  data  into  the  interior  of  the  model,  we  conducted  a  sim¬ 
ple  predictability  experiment.  Travel  time  data  was  withheld  from  the  interior, 
and  only  the  boundary  of  the  model  was  updated.  Here,  we  wished  to  demon¬ 
strate  how  long  a  time  period  the  model  was  able  to  continue  to  predict  the 
interior  travel  time  data.  Further,  we  wished  to  estimate  how  often  interior 
data  needed  to  be  assimilated  into  the  model.  To  evaluate  at  what  point  the 
model  diverges  from  the  measurements,  we  compared  the  differences  between 
data  that  the  boundary  driven  model  predicts  and  the  AMODE  data.  Panel  A 
in  Figure  4.6  shows  the  RMS  level  of  the  these  residuals  and  the  RMS  level  of 
the  AMODE  travel  times.  After  a  period  of  about  60  days,  the  predicted  data 
has  diverged  from  the  measurements.  This  is  an  indication  that  the  QG  model 
requires  interior  updating  within  at  least  the  60  day  time  frame  in  order  to 
maintain  any  degree  of  predictability.  The  wavenumber  spectra  of  the  model 
are  shown  in  Panels  C  and  D  of  Figure  4.6.  They  show  the  energy  to  be  pre¬ 
dominantly  low  wavenumber,  with  the  barotropic  mode  having  somewhat  more 
temporal  variability. 

Results  of  the  interior  assimilation  are  shown  in  Figure  4.7.  The  x2  of  the  fit 
grows  a  bit  during  the  220  day  assimilation  period,  but  remains  within  reason. 
The  RMS  of  the  innovations  consistently  remains  lower  than  than  the  RMS  of 
the  data.  Panel  B  shows  the  final  baroclinic  assimilation  map  for  year  day  350. 

Panels  C  and  D  show  the  wavenumber  spectra  for  the  model.  The  wavenum¬ 
ber  spectra  of  both  modes  are  red.  Comparing  with  Figures  4.3  and  4.5,  the 
wavenumber  estimates  are  slightly  different  than  those  estimated  using  the 
Rossby  wave  models  which  tended  to  have  more  energy  in  higher  wavenumber 
components. 

Plots  of  the  time  evolution  of  the  stream  function  and  vorticity  RMS  ener¬ 
gies  and  the  kinetic  energies  are  shown  in  Figure  4.8.  All  of  the  energy  levels 
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are  initially  low,  and  begin  to  increase  after  year  day  180.  This  general  in¬ 
crease  in  energy  also  matches  an  increase  in  the  baroclinic  energy  estimated  in 
El  using  the  Rossby  wave  model. 

Panel  A  shows  the  stream  function  RMS  energy  as  a  function  of  time  for 
each  of  the  six  levels  in  the  QG  model  (  “1”  being  shallowest  to  “6”  the  deepest). 
The  upper  two  levels  contain  most  of  the  energy,  roughly  twice  that  of  the  lower 
4  levels.  Panel  B  shows  the  RMS  energies  of  the  stream  function  in  terms  of 
modes.  Nearly  all  of  the  variability  is  contained  in  the  first  two  modes.  Some 
energy  is  being  transfered  to  higher  modes  by  the  model,  as  can  be  seen  by  the 
increase  in  the  energy  for  the  third  mode.  This  transfer  of  energy  is  a  non-linear 
process,  and  arises  from  mode-mode  interactions  within  the  QG  model.  Panels 
C  and  D  show  level/mode  time  histories  for  the  vorticity.  They  show  a  similar 
story  as  the  stream  function:  most  of  the  energy  is  contained  in  the  upper  levels 
and  first  two  modes,  but  a  bit  more  energy  appears  to  be  imparted  to  the  higher 
vorticity  modes  (Panel  D).  The  upper  panel  shows  the  RMS  kinetic  energy  in 
terms  of  both  level  and  modes.  The  upper  level  kinetic  energies  appear  to  be 
between  20-30  (cm2  s~2),  and  the  kinetic  energy  is  roughly  equally  partitioned 
between  the  first  two  modes. 
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Figure  4.1:  Panel  A)  RMS  error  levels.  Panel  B)  Final  assimilation  and  true 
(control)  field  estimates  for  the  baroclinic  mode.  Panel  C)  Differences  be¬ 
tween  the  estimate  and  control  ocean  for  the  barotropic  mode  in  terms  of  the 
wavenumber  space  (CpMm).  D)  Differences  between  the  estimate  and  control 
ocean  for  the  baroclinic  mode  in  terms  of  the  wavenumber  space.  Panel  E)  Es¬ 
timated  spectral  errors  in  the  barotropic  mode.  Panel  F)  Estimated  spectral 
errors  in  the  baroclinic  mode. 
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Figure  4.2:  Comparison  between  error  in  the  wavenumber  estimate  as  a  func¬ 
tion  of  time  for  the  point  measurements  (dashed  lines)  and  the  tomographic 
measurements  (solid  lines).  This  figure  shows  the  error  reduction  in  wavenum¬ 
ber  space  at  3  different  times  during  the  assimilation  of  the  point  and  tomo¬ 
graphic  data. 
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Figure  4.3:  Panel  A  )  RMS  error  levels.  Dark  yellow/blue  are  the  residuals 
of  the  forecast/analysis.  Panel  B)  Final  assimilation  for  the  baroclinic  mode. 
Panel  C)  Barotropic  mode  amplitude  in  wavenumber  space  (CpMm).  D)  Baro¬ 
clinic  mode  amplitude  in  wavenumber  space.  Panel  E)  Estimated  spectral  er¬ 
rors  in  the  barotropic  mode.  Panel  F)  Estimated  spectral  errors  in  the  baroclinic 
mode. 
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Figure  4.4:  Panel  A)  RMS  error  levels.  Dark  yellow/blue  are  the  residuals  of 
the  forecast/analysis,  light  yellow/blue  are  the  errors  of  the  magnitudes  of  the 
bc/bt  velocity  estimate.  Panel  B)  Final  assimilation  estimates  for  the  baroclinic 
mode  in  physical  space.  Panel  C)  Differences  between  the  estimate  and  control 
ocean  for  the  barotropic  mode  in  terms  of  the  wavenumber  space  (CpMm).  D) 
Differences  between  the  estimate  and  control  ocean  for  the  baroclinic  mode  in 
terms  of  the  wavenumber  space.  Panel  E)  Estimated  spectral  errors  in  the 
barotropic  mode.  Panel  F)  Estimated  spectral  errors  in  the  baroclinic  mode. 
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Figure  4.5:  Panel  A)  RMS  error  levels. 

Dark  yellow/blue  are  the  residuals  of  the  forecast/analysis,  light 
yellow/blue  are  the  errors  of  the  magnitudes  of  the  bc/bt  velocity  estimate. 
Panel  B)  Estimated  advection  velocities  and  error  bars  (shaded). 

Panel  C)  Differences  between  the  estimate  and  control  ocean  for  the  barotropic 
mode  in  terms  of  the  wavenumber  space  (CpMm).  Panel  C)  Barotropic  mode 
amplitude  in  wavenumber  space.  D)  Baroclinic  mode  amplitude  in  wavenum¬ 
ber  space.  Panel  E)  Estimated  spectral  errors  in  the  barotropic  mode.  Panel  F) 
Estimated  spectral  errors  in  the  baroclinic  mode. 
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Figure  4.6:  Panel  A)  RMS  of  the  measured  data  and  the  difference  with  data 
predicted  by  the  boundary  force  HOOM  model.  Panel  B)  Final  field  in  physical 
space.  Panel  C)  Amplitude  of  the  barotropic  mode  in  spectral  domain  (CpMm). 
Panel  D)Amplitude  of  the  baroclinic  mode  in  spectral  domain. 
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Figure  4.7:  Panel  A)  RMS  errors.  Panel  B)  final  analysis  field.  Panel  C) 
Barotropic  mode  amplitude  in  wavenumber  space.  D)  Baroclinic  mode  ampli¬ 
tude  in  wavenumber  space  (CpMm).  Panel  E)  Estimated  spectral  errors  in  the 
barotropic  mode.  Panel  F)  Estimated  spectral  errors  in  the  baroclinic  mode. 
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Figure  4.8:  Quasi-geostrophic  energies.  (See  text  for  description.) 
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Chapter  5 

DISCUSSION 

We  began  with  the  hypothesis  that  acoustic  tomographic  measurements  are 
uniquely  important  because  they  are  sensitive  to  low  wavenumber  properties  of 
the  ocean,  and  that  the  inherent  integrating  nature  of  the  measurements  must 
be  accounted  for  when  the  data  are  assimilated  into  numerical  ocean  models. 

Munk  and  Wunsch  (1979)  first  introduced  the  concept  of  ocean  tomography 
as  a  means  of  determining  the  eddy-scale  structure  of  the  ocean.  They  showed 
that  acoustic  travel  times  could  be  used  to  measure  integral  or  average  proper¬ 
ties  of  soundspeed  (which  to  first  order  is  linearly  related  to  temperature)  and 
current  velocity  along  raypaths  between  an  array  of  sources  and  receivers.  Fur¬ 
ther,  they  showed  to  first  order  that  the  one-way  travel  times  could  be  related 
to  the  integrated  temperature  anomaly  along  the  ray,  and  that  the  differenced 
travel  times  could  be  related  to  the  integrated  current  velocity  along  the  ray. 

By  using  an  array  of  acoustical  transceivers  these  temperature  and  current 
fields  can  be  mapped  over  a  region  of  the  ocean.  Although  the  “times  of  flight” 
of  the  acoustical  signals  between  transceivers  are  nearly  instantaneous,  typi¬ 
cal  tomographic  experiments  are  designed  to  repeatedly  sample  the  ocean  and 
generate  long  time  series  of  measurements  over  several  weeks  or  months.  Dur¬ 
ing  these  time  periods  the  ocean  fields  may  appreciably  change  and  the  latency 
or  time  dependence  of  the  data  must  be  taken  into  account. 

Over  the  last  two  decades,  a  growing  number  of  tomographic  experiments 
have  been  conducted  in  many  of  the  world’s  oceans.  One  such  experiment, 
AMODE,  was  located  north  of  Puerto  Rico.  During  this  experiment,  an  ex¬ 
traordinarily  long  time  series  (220  days)  of  tomographic  data  was  collected. 
These  data  were  specifically  designed  to  resolve  the  eddy  scale  temperature 
and  current  variability. 

In  part,  this  research  was  motivated  by  the  availability  and  need  to  ana¬ 
lyze  a  data  set  such  as  the  AMODE  data,  as  well  as  by  a  series  of  workshops 
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sponsored  by  ONR.  The  primary  goal  of  these  workshops  was  to  foster  com¬ 
munication  between  the  ocean  modeling  community  and  the  data  acquisition 
community.  Naturally,  integral  measurements  such  as  tomography  are  only 
one  piece  of  a  total  ocean  observing  system.  These  workshops  discussed  how 
to  incorporate  many  different  forms  of  data  into  time-dependent  ocean  models 
including  point  or  localized  measurements. 

5.1  Parameterization 

In  order  to  explain  why  we  parameterized  our  model  the  way  we  did,  we  need 
to  step  back  and  look  at  the  fundamental  properties  of  the  tomography  data 
which  measure  integral  properties  along  acoustical  ray  paths.  First,  we  con¬ 
sidered  a  set  of  raypaths  along  any  given  vertical  slice  between  a  source  and 
receiver.  These  raypaths  are  determined  by  the  soundspeed  structure  in  the 
ocean  which  acts  as  a  dispersive  wave  guide.  For  purposes  of  this  research 
we  characterized  these  raypaths  using  an  average  or  background  soundspeed 
field  and  we  estimated  perturbations  to  this  reference.  For  typical  soundspeed 
anomalies  associated  with  the  eddy  fields  in  the  AMODE  region,  the  acous¬ 
tical  paths  were  assumed  fixed,  i.e.,  their  locations  were  unperturbed  by  the 
soundspeed  anomaly.  This  assumption  was  tested  by  calculating  ray  trajecto¬ 
ries  through  perturbations  that  were  reasonable  for  this  region.  This  assump¬ 
tion  was  found  to  be  valid. 

Because  these  rays  spend  most  of  the  time  near  their  turning  depths, 
they  are  particularly  sensitive  to  temperature  and  current  anomalies  at  those 
depths.  The  raypaths  are  controlled  by  the  dispersive  wave-guide  properties 
of  the  ocean  propagation.  Only  a  limited  number  are  typically  discernible 
between  the  source  and  receiver  (15-20  for  the  AMODE  experiment).  Hence 
along  any  vertical  tomographic  slice,  the  rays  will  only  sample  a  limited  num¬ 
ber  of  depths  and  produce  rather  coarse  vertical  resolution.  In  the  horizontal 
plane,  resolution  is  determined  by  the  number  of  crossing  paths  between  the 
sources  and  receivers.  Theoretically,  we  could  have  an  unlimited  number  of 
them,  but  because  of  limited  budgets,  we  must  make  do  with  as  few  as  possi¬ 
ble. 

Given  these  limitations  of  the  tomographic  vertical  and  horizontal  sampling 
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properties,  some  method  had  to  be  developed  to  parameterize  soundspeed  and 
current  fields  in  the  horizontal  and  vertical  directions.  The  parameterization 
had  to  serve  three  purposes:  1)  reduce  the  number  of  model  parameters  which 
we  needed  to  fit  the  travel  time  data,  2)  provide  a  means  of  interpolation,  and  3) 
provide  a  means  of  converting  between  soundspeed  and  current  parameteriza¬ 
tion  and  the  parameterization  used  by  the  circulation  models  (stream  function 
and  vorticity  which  will  be  discussed  shortly). 

There  are  several  possible  basis  sets  which  could  be  used.  For  instance, 
triangle  or  Gaussian  shaped  functions  (these  were  used  early  on  in  this  re¬ 
search,  but  discarded  for  reasons  stated  below)  or  functions  that  are  based 
on  the  statistical  soundspeed/current  variability  of  the  circulation  model  or 
the  ocean  itself  (EOF’s).  Because  EOFs  are  eigenvectors  of  the  fields  covari¬ 
ance  matrix,  they  provide  the  most  efficient  representation  of  the  field  with  the 
fewest  modes.  Unfortunately,  EOF’s  do  not  have  continuous  derivatives  which 
is  a  necessary  condition  to  produce  realistic  soundspeed  modes.  In  other  words, 
rays  could  not  be  re-traced  through  the  soundspeed  field.  In  this  research,  we 
ultimately  chose  a  parameterization  related  to  the  physics  of  the  problem. 

In  the  vertical,  we  chose  to  use  dynamical  modes.  These  modes  are  based 
on  QG  theory,  and  only  the  first  two  modes  are  necessary  to  explain  98%  of  the 
variance  in  the  temperature  and  current  fields  of  the  AMODE  region.  Not  only 
are  they  related  to  the  physics  of  the  flow,  but  they  are  an  efficient  basis  set  as 
well.  Also,  these  modes  can  be  conveniently  related  to  the  soundspeed  modes 
using  buoyancy  frequency  and  the  gradient  of  potential  temperature. 

With  the  barotropic  mode  having  no  zero  crossings,  the  first  baroclinic  mode 
has  a  crossing  at  a  depth  of  about  1400  m  in  the  AMODE  region.  The  barotropic 
mode  is  constant  in  depth  and  has  no  temperature  or  soundspeed  mode  signa¬ 
ture,  but  it  does  have  a  significant  current  velocity  signature.  This  mode  was 
used  to  describe  the  vertically  uniform  current  field.  Theoretically,  in  the  ab¬ 
sence  of  bottom  topography,  its  radius,  or  length  scale  of  deformation,  is  infi¬ 
nite. 

The  first  baroclinic  mode  has  a  single  zero  crossing  and  a  maximum  at  the 
surface.  Its  equivalent  soundspeed  mode  has  a  maximum  at  a  depth  of  about 
700  m  and  a  theoretical  deformation  radius  of  about  50  km.  Each  successive 
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mode  will  have  an  additional  zero  crossing  and  a  shorter  deformation  radius. 

In  the  horizontal,  we  chose  to  use  sines  and  cosines.  This  choice  was  based 
on  two  main  factors.  First,  we  wished  to  estimate  the  spectrum  of  Rossby 
waves.  Secondly,  based  on  our  hypothesis,  we  wished  to  evaluate  the  per¬ 
formance  of  the  data  assimilation  experiments  in  terms  of  the  wavenumber 
domain.  We  previously  mentioned  that  each  vertical  dynamical  mode  has  a 
corresponding  radius  of  deformation;  it  is  important  to  keep  in  mind  that  the 
model  wavenumber  domain  should  be  made  large  enough  to  resolve  the  partic¬ 
ular  mode  of  interest.  If  higher  modes  are  used,  more  wavenumbers  should  be 
included  in  the  wave  field  in  order  to  resolve  the  higher  order  modes. 

A  drawback  to  using  sines  and  cosines  is  that  they  do  a  very  poor  job  of  ex¬ 
trapolation  away  from  the  data  domain.  Also,  they  can  force  periodic  boundary 
conditions  on  the  solution  if  the  physical  domain  boundaries  are  not  at  least 
1.5-2  times  larger  than  the  data  domain.  In  this  research,  we  strived  to  elimi¬ 
nate  this  effect  by  choosing  the  appropriate  model  domains. 

We  do  in  effect  extrapolate  the  Rossby  wave  field  in  order  to  produce  bound¬ 
ary  conditions  for  the  QG  model;  however,  we  try  to  minimize  the  extrapolation 
errors  on  the  interior  domain  of  the  model  by  choosing  the  QG  domain  to  be 
sufficiently  large  so  that  the  assimilation  of  interior  data  can  correct  this. 

5.2  Twin  Experiments 

We  made  use  of  several  twin  experiments.  The  utility  of  the  twin  experiment 
was  to  demonstrate  the  effectiveness  of  a  particular  data  set  or  combination 
of  data  sets  when  they  are  assimilated  or  combined  with  an  ocean  circulation 
model.  The  main  purpose  of  the  twin  experiments  was  to  make  comparisons 
between  a  control  model  and  another  model  which  is  being  updated  or  corrected 
using  synthetic  data  derived  from  the  control  ocean.  Data  assimilation  was 
used  to  make  the  corrections  to  the  twin  model’s  parameters.  The  twin  model 
was  advanced  in  time  using  the  same  physics  as  the  control  model  (this  does 
not  imply  that  the  model  was  noise  free  because  the  interior  and  boundary 
are  not  exact).  Hypothetically,  if  the  data  were  perfect  (noise  free  and  able  to 
resolve  all  scales  of  the  model),  the  model  noise  estimates  were  perfect,  and  an 
assimilation  method  that  fully  propagated  the  estimates  and  their  errors  were 
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used,  results  of  the  two  models  should  converge.  To  the  degree  that  any  of  these 
things  is  less  than  perfect,  the  assimilation  and  the  truth  will  not  converge. 

Typically,  there  is  little  point  in  conducting  a  twin  experiment  that  uses 
perfect  data,  so  we  sampled  the  control  model  and  contaminated  the  synthetic 
observations  with  noise  commensurate  with  some  a  priori  noise  assumptions 
( R ).  The  twin  experiment  was  used  to  determine  how  to  adjust  Q  (the  model 
forecast  covariance)  so  that  the  estimate  and  the  control  oceans  converge  opti¬ 
mally.  In  as  much  as  we  have  confidence  that  this  twin  experiment  represented 
nature,  we  used  these  estimates  for  Q  to  guide  us  when  assimilating  real  data. 

Although  we  learned  a  great  deal  from  the  twin  experiments,  for  example, 
sensitivity  of  the  estimate  to  Q,  there  are  several  things  we  could  not  learn. 
For  example,  they  cannot  tell  us  if  our  measurement  model,  h(x),  is  correctly 
modeling  our  observations.  The  measurements  will  always  be  self  consistent 
with  the  model.  The  normalized  x2  estimate  (the  distribution  of  the  residuals 
weighted  by  their  errors)  from  the  fit  will  be  sensitive  to  changes  in  Q,  but  not 
to  changes  in  the  data  error,  since  the  noise  added  to  the  data  is  statistically 
consistent  with  the  a  priori  data  error.  This  means  that  in  order  to  validate  the 
measurement  model,  separate  comparisons  need  to  be  made  between  modeled 
data  and  actual  data. 

5.3  Assimilation  Techniques 

Data  assimilation  is  a  term  used  for  various  techniques  of  incorporating  or 
combining  data  into  a  time  dependent  model.  Intrinsically,  all  of  the  assimi¬ 
lation  methods  can  be  considered  simplifications  or  recursive  estimates  of  the 
more  general  least  squares  estimate  of  the  model  state  over  all  space  and  time. 
Except  for  all  but  the  simplest  linear  models  (even  with  the  parameterization 
mentioned  previously),  this  quickly  becomes  computationally  infeasible.  There¬ 
fore  sequential  methods,  such  as  the  Kalman  filter,  are  used  to  estimate  the 
model  state  and  its  errors  at  successive  time  steps  during  the  evolution  of  the 
model. 

The  Kalman  filter  requires  in  its  transition  step  the  computation  of  the 
mode  error  covariance,  which  is  very  computationally  intensive.  For  linear 
models  (the  Rossby  wave  models),  the  transition  matrix  can  be  thought  of  as 
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a  tri-diagonal  matrix  which  phase  shifts  the  wavenumber  components  of  the 
model.  The  full  Kalman  filter  was  used  for  the  Rossby  wave  models. 

However,  for  the  QG  model,  there  is  no  simple  expression  for  the  transition 
matrix.  The  transition  matrix  is  dependent  upon  the  current  estimate  of  the 
stream  function  and  vorticity  fields.  For  this  model,  optimal  interpolation  (01) 
was  used,  wherein  the  model  was  used  to  forecast  the  state  at  the  appropriate 
time  step,  but  not  its  error  covariance.  This  can  be  viewed  as  using  I,  the  iden¬ 
tity  matrix,  in  place  of  the  transition  matrix  for  the  calculation  of  the  forecast 
model  error  covariance.  Although  advection  of  information  is  lost  in  such  a  co- 
variance  update,  information  about  correlations  in  the  data  is  retained  (from 
the  analysis  step  Pa  =  (I  -  KH)Pf). 

A  number  of  twin  experiments  were  conducted  using  the  linear  Rossby  wave 
model  in  order  to  evaluate  the  penalty  associated  with  this  sub-optimal  tran¬ 
sition  step.  After  tomographic  data  was  assimilated  over  several  realizations 
of  the  wave  field,  the  RMS  error  in  the  Kalman  filter  estimate  was  10%,  and 
the  RMS  error  in  the  01  estimate  was  18%.  The  difference  between  the  two 
eventually  became  negligible  after  assimilating  more  data. 

5.4  Errors 

The  twin  experiments  used  theoretical  travel  times  contaminated  with  random 
errors  whose  variances  were  simply  specified  as  constants  for  the  sum  travel 
times  and  the  difference  travel  times.  These  constants  were  based  on  average 
error  estimates  obtained  from  experience  with  the  AMODE  data.  The  constant 
used  for  the  reciprocal  travel  times,  a  =  ±3  ms,  was  larger  than  that  normally 
observed  for  experiments  such  as  this  (Dushaw  et  al.,  1996).  It  is  believed 
that  these  errors  may  be  reduced  by  improved  mooring  tracking  analysis  of  the 
AMODE  data  (B.  Dushaw,  personal  communication,  May  1999). 

The  error  estimates  contained  the  accumulative  effects  of:  1)  residual  moor¬ 
ing  motion  (sum  travel  times  only),  2)  uncorrected  clock  drift,  and  3)  internal 
wave  noise.  Additional  errors  from  assuming  that  the  raypaths  were  unper¬ 
turbed  are  on  the  order  0.1  ms  and  were  not  included  in  the  error  budget.  In 
hindsight  this  error  source  should  have  been  added  for  completeness.  In  the 
twin  simulations,  the  travel  time  errors  were  artificially  lowered  (by  dividing 
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the  error  by  \/5)  in  an  effort  to  simulate  the  use  of  15  ray  arrivals  while  using 
only  3.  This  was  valid  given  that  we  were  only  using  two  vertical  modes,  which 
were  adequately  constrained  (the  eigenvalue  spectrum  only  contains  2  signif¬ 
icant  eigenvalues).  Therefore,  if  the  experiments  were  repeated  using  all  15 
available  rays,  the  results  should  not  be  significantly  different. 

5.5  Models  and  Results 

We  chose  to  assimilate  the  tomographic  data  using  three  different  models.  The 
two  Rossby  wave  models  were  specifically  chosen  for  their  simple  dynamics. 
These  simple  dynamics  allowed  the  Kalman  filter  to  be  used  for  the  updating. 
We  did  not  expect  them  to  realistically  describe  the  flow  in  the  AMODE  region. 
The  QG  model  was  chosen  because  it  was  a  non-linear  extension  of  the  linear 
Rossby  wave  models  (the  Jacobian  and  bottom  topography  set  them  apart)  and 
provided  a  more  realistic  situation. 

The  results  of  the  twin  experiments  with  the  linear  Rossby  models  showed 
that  the  tomographic  data  (El)  resolved  the  low  wavenumber  components  of 
the  Rossby  wave  field.  The  RMS  difference  between  the  final  assimilation  es¬ 
timate  of  the  Rossby  wave  field  and  the  control  ocean  at  the  end  of  the  assimi¬ 
lation  period  was  decreased  by  about  50%  from  its  starting  value.  We  explored 
using  various  estimates  for  Q.  These  included  using  an  adaptive  Q,  and  <5  =  0 
at  all.  Although  x2  was  sensitive  to  these  changes,  the  final  RMS  difference  was 
not.  Maps  of  the  difference  between  the  fields  did  show  some  minor  improve¬ 
ments  when  Q  =  0.  This  is  expected  because  with  non-zero  Q  we  effectively 
over-weighted  the  observations  and  under-weighted  the  model. 

Point  data  were  also  assimilated  into  the  linear  Rossby  wave  model.  The 
number  of  measurements  was  approximately  equal  to  the  number  of  sum  travel 
times.  The  data  error  for  the  point  measurements  was  adjusted  so  that  the  av¬ 
erage  error  in  the  final  wavenumber  estimates  at  the  end  of  the  assimilation 
period  was  comparable  to  that  produced  by  the  tomographic  assimilation.  Al¬ 
though  the  error  wavenumber  spectra  estimates  for  both  data  were  initially 
similar,  at  the  end  of  the  300  day  period  the  tomographic  errors  for  the  low 
wavenumber  components  were  less  than  those  for  the  point  measurements. 
This  finding  was  consistent  with  the  integral  nature  of  the  tomographic  mea- 
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surements.  The  noise  level  of  the  point  measurements  was  effectively  ±90  ms, 
about  20  times  larger  than  the  travel  time  errors. 

Results  of  the  twin  experiment  that  included  advection  were  nearly  iden¬ 
tical  to  the  non-advection  experiment.  The  RMS  reduction  in  the  difference 
between  the  estimated  and  true  stream  function  fields  decreased  by  50%. 
The  barotropic  advection  components,  which  are  determined  by  the  difference 
travel  times,  could  not  be  resolved.  Most  likely  this  is  because  the  very  small 
travel  time  signature  of  this  barotropic  advective  component  («  0.8  ms)  was 
swamped  by  the  difference  travel  time  noise. 

5.6  Real  Data 

Both  of  the  Rossby  wave  models  fit  the  actual  AMODE  travel  time  data  in 
a  manner  consistent  with  the  assumed  error  levels  in  the  data  and  a  priori 
model  errors.  The  correlation  coefficient  between  the  measured  data  and  that 
predicted  by  the  models  was  0.84,  and  the  variance  reduction  was  61%  for  both 
models,  indicating  that  both  models  on  average  fit  the  data  equally  as  well. 

The  normalized  y2  estimates  from  both  model  fits  were  nearly  identical,  and 
remained  near  1  until  year  day  300,  when  they  increased.  This  corresponds  to 
the  time  at  which  most  of  the  transceivers  failed  and  there  was  only  data  from 
2  or  3  moorings.  Since  the  wavefields  estimated  by  both  models  were  nearly 
identical,  we  describe  next  the  wavenumber  fit  from  the  non-advective  model. 

The  spectra  of  both  modes  had  roughly  the  same  shape.  This  is  not  un¬ 
expected,  since  the  same  a  priori  spectrum  was  assumed  for  both  modes.  In 
hindsight,  this  was  not  wise.  The  barotropic  spectrum  should  have  been  redder 
(based  on  the  internal  deformation  radii). 

Some  of  the  details  of  the  time  evolution  of  the  mode  spectra  are  quite  in¬ 
teresting.  For  instance,  there  were  several  episodic  rises  and  declines  in  the 
strengths  of  the  two  modes.  These  rises  and  declines  did  not  occur  at  the  same 
time  for  both  modes.  It  is  suggested  that  these  packets  of  activity  represent 
periods  of  more  intense  wave  activity.  For  the  barotropic  waves,  there  were 
two  brief  increases  near  year  days  100  and  160,  and  then  a  longer  duration  in¬ 
crease  between  year  days  250  and  350.  For  the  baroclinic  mode,  there  was  a  60 
day  intensification  starting  near  year  day  180  which  was  followed  by  a  slight 
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decrease  in  amplitude. 

The  fits  to  the  QG  model  reduced  somewhat  more  of  the  travel  time  variance 
(66%)  than  the  Rossby  wave  models  (61%).  The  x2  of  the  fit  was  similar  to  the 
Rossby  wave  fits,  with  a  similar  increase  after  day  300.  The  RMS  level  of  the 
forecast  residuals  was  greater  than  for  the  Rossby  wave  model.  We  feel  we 
could  improve  this  with  better  estimates  of  Q. 

The  time  evolution  of  the  mode  spectrums  was  a  bit  different  from  the 
Rossby  wave  fits.  The  barotropic  mode  showed  an  increase  in  amplitude  near 
year  day  200,  about  50  days  earlier.  Also  this  intense  period  only  lasted  about 
50  days.  Another  intense  period  occurred  near  year  day  280  and  died  off  by  day 
300.  The  evolution  of  the  baroclinic  mode  spectra  shows  no  major  difference 
with  the  Rossby  wave  evolution.  However,  examination  of  the  physical  space 
maps  of  the  two  did  reveal  subtle  differences  in  the  amplitude  and  phase  of  the 
stream  function  fields. 

The  spectrum  of  the  baroclinic  mode  for  the  QG  models,  on  average,  has 
slightly  more  high  wavenumber  energy  than  the  spectra  of  the  linear  mod¬ 
els.  This  higher  wavenumber  energy  is  a  result  of  the  model’s  interaction  with 
topography.  More  energy  in  these  higher  wavenumbers  may  account  for  the 
increased  variance  reduction  in  the  QG  model’s  fit. 

Regional  models,  such  as  HOOM  (HOOM  is  the  Harvard  Open  Ocean  QG 
model  used  in  this  research),  require  boundary  conditions,  and  some  means  of 
extrapolation  must  be  used  if  the  interior  data  is  insufficient  (in  time  or  space) 
to  supply  those  conditions  (Bennett,  1992).  There  are  a  number  of  possible  al¬ 
ternatives  besides  the  method  that  was  used  here  (Kalman  filtered  estimates 
of  Rossby  waves),  such  as  extrapolation  from  the  interior  data  or  large-scale 
forcing  from  global  circulation  models.  This  is  conceptually  attractive,  but  as 
data  is  assimilated  into  the  the  model’s  interior,  serious  discontinuities  could 
develop  near  the  boundary  region.  Secondly,  the  boundaries  could  have  been 
expanded  outward  several  hundred  kilometers  from  the  assimilation  region. 
This  is  numerically  unattractive  because  it  requires  extra  computation  time.  A 
third  possibility  for  supplying  boundary  conditions  would  be  to  use  data  from 
the  interior  in  conjunction  with  some  form  of  objective  mapping.  This  was  at¬ 
tempted  early  in  this  research  but  found  to  be  problematic  because  of  sudden 
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temporal  jumps  in  the  boundary,  which  the  circulation  model  did  not  have  time 
to  adjust  for.  These  jumps  resulted  in  spurious  growth  in  the  stream  function 
and  vorticity  fields. 

The  approach  to  boundary  conditions  used  in  the  QG  model  experiment  at¬ 
tempted  to  make  use  of  the  interior  data  and  still  prescribe  the  boundary  condi¬ 
tions  in  a  manner  consistent  with  the  model  dynamics  and  the  interior  data.  In 
this  research,  we  re-fitted  the  interior  data  once;  this  process  could  be  further 
improved  with  more  iterations.  That  is,  after  the  first  iteration,  we  could  have 
refined  the  boundary  conditions  and  re-assimilated  the  interior  data  a  second 
time. 

Experiment  E3  indicated  that  the  QG  model  fit  more  data  variance  than 
either  of  the  Rossby  wave  models.  Fundamentally,  this  was  because  the  QG 
model  offers  more  degrees  of  freedom  than  does  the  Rossby  wave  model.  In 
essence  the  QG  model  allowed  for  a  “richer”  spectrum  of  waves  in  its  interior 
flow  field  compared  to  the  limited  number  of  spectral  coefficients  offered  by  the 
Rossby  wave  model.  Non-linear  effects  of  the  QG  model  may  also  have  been  a 
factor.  The  non-linear  transfer  of  energy  to  the  third  baroclinic  mode  shown  in 
Experiment  E3  may  have  accounted  for  the  QG  model  capturing  more  of  the 
data  variance. 

E3  also  indicated  important  differences  between  the  time  evolution  of  the 
waves  that  were  fit  using  the  Rossby  wave  model  and  the  QG  model.  Although 
the  spectrum  of  the  baroclinic  mode  appeared  similar  in  both  model  estimates, 
the  barotropic  mode  had  more  temporal  variability  in  the  QG  esimate.  This 
may  be  attributed  to  the  additional  topography  and  the  non-linear  effects  in 
the  QG  model.  It  would  be  useful  in  the  future  to  examine  these  differences 
between  the  linear  and  non-linear  fits  more  fully. 

The  RMS  of  the  innovation  sequence  was  considerably  higher  than  the  RMS 
of  the  residuals.  This  was  likely  due  to  mis-specifying  the  model  errors,  Q.  As¬ 
sumptions  about  the  data  and  model  errors  are  important  factors  in  reducing 
the  innovation  sequence.  For  instance,  underestimating  Q  could  cause  the  fore¬ 
cast  covariance  to  collapse  and  result  in  the  measurements  being  ignored  in  the 
assimilation.  Conversely,  overestimating  Q  would  cause  the  measurements  to 
be  weighted  too  heavily.  A  more  preferable  estimate  for  Q  could  be  obtained 
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from  a  twin  experiment  run  using  the  QG  model  itself. 

This  work  is  just  the  beginning  of  what  could  be  accomplished  with  the 
AMODE  data  set.  We  have  a  stronger  understanding  of  how  to  incorporate 
tomographic  data  into  circulation  models  and  the  types  of  results  one  can  ex¬ 
pect.  Further,  we  have  gained  insight  into  the  usefulness  of  assimilating  tomo¬ 
graphic  measurements  into  numerical  models.  There  is  great  opportunity  to 
extend  this  research  to  learn  even  more  about  the  AMODE  region  and,  at  the 
same  time,  learn  more  about  the  subtleties  of  assimilating  tomographic  data 
into  numerical  models. 
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